Rice-Wheat System Spatial Decision Support Tool

Author

Maxwell Mkondiwa (m.mkondiwa@cgiar.org)

1 Introduction

We present a rice-wheat systen spatial decisionb support tool which relies on data from the matched landscape crop assessment surveys in eastern India and a multivariate geoadditive Bayesian model to predict entry points for system optimization in the area of interest.

2 Agronomic Functional Relationships: Crop Specific Regressions

2.1 Rice

Code 
library(sp)
library(mgcv)
library(bamlss)

# Multivariate geoadditive model
# remotes::install_git("https://git.uibk.ac.at/c4031039/mvnchol")
# 
library(mvnchol)
library(BayesX)
library(R2BayesX)
library(sf)
library(spdep)
library(rio)
Irrig_Rev_rice_wheat <- import("data/Irrig_Rev_rice_wheat_Updated.csv")

Irrig_Rev_rice_wheat <- subset (Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$Longitude)))
Irrig_Rev_rice_wheat <- subset (Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$Latitude)))

library(janitor)
# Irrig_Rev_rice_wheat <- clean_names(Irrig_Rev_rice_wheat)
library(tidyr)
#Irrig_Rev_rice_wheat <- Irrig_Rev_rice_wheat %>% drop_na()

# library(lubridate)
# library(anytime)
# Irrig_Rev_rice_wheat$januaryfirst2017 <- ymd("2017-01-01")

# Irrig_Rev_rice_wheat$sowdate_fmt_wheat_day <- Irrig_Rev_rice_wheat$sowdate_fmt - Irrig_Rev_rice_wheat$januaryfirst2017

# Irrig_Rev_rice_wheat$sowdate_fmt_rice_day <- Irrig_Rev_rice_wheat$sowdate_fmt_rice - Irrig_Rev_rice_wheat$januaryfirst2017

Irrig_Rev_rice_wheat$harvest_day_rice <- Irrig_Rev_rice_wheat$l_crop_duration_days_rice + Irrig_Rev_rice_wheat$sowdate_fmt_rice_day

# Irrig_Rev_rice_wheat$harvest_day_wheat <- Irrig_Rev_rice_wheat$l_crop_duration_days + Irrig_Rev_rice_wheat$sowdate_fmt_wheat_day


shpname <- file.path(getwd(), "shp", "India_aoi_sf_sp")

India_aoi_sp_bnd <- BayesX::shp2bnd(shpname = shpname, regionnames = "District", check.is.in = F)
Reading map ... finished
Note: map consists originally of 50 polygons
Note: map consists of 47 regions
Code 
f_rice_yield_MRF <- list(
    b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + s(sowdate_fmt_rice_day) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    sigma ~ 1 + rice_duration_class_long + s(sowdate_fmt_rice_day) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re")
)



K <- neighbormatrix(India_aoi_sp_bnd)
head(K)
           Araria Arwal Aurangabad Banka Begusarai Bhagalpur Arah Buxar
Araria          4     0          0     0         0         0    0     0
Arwal           0     6         -1     0         0         0   -1     0
Aurangabad      0    -1          3     0         0         0    0     0
Banka           0     0          0     3         0        -1    0     0
Begusarai       0     0          0     0         5         0    0     0
Bhagalpur       0     0          0    -1         0         6    0     0
           Darbhanga Gaya Gopalganj Jamui Jehanabad Kaimur Katihar Khagaria
Araria             0    0         0     0         0      0       0        0
Arwal              0   -1         0     0        -1      0       0        0
Aurangabad         0   -1         0     0         0      0       0        0
Banka              0    0         0    -1         0      0       0        0
Begusarai          0    0         0     0         0      0       0       -1
Bhagalpur          0    0         0     0         0      0      -1       -1
           Kishanganj Lakhisarai Madhepura Madhubani Munger Muzaffarpur Nalanda
Araria             -1          0        -1         0      0           0       0
Arwal               0          0         0         0      0           0       0
Aurangabad          0          0         0         0      0           0       0
Banka               0          0         0         0     -1           0       0
Begusarai           0         -1         0         0     -1           0       0
Bhagalpur           0          0        -1         0     -1           0       0
           Nawada WestChamparan Patna EastChamparan Purnia Rohtas Saharsa
Araria          0             0     0             0     -1      0       0
Arwal           0             0    -1             0      0     -1       0
Aurangabad      0             0     0             0      0     -1       0
Banka           0             0     0             0      0      0       0
Begusarai       0             0    -1             0      0      0       0
Bhagalpur       0             0     0             0     -1      0       0
           Samastipur Saran Sheikhpura Sheohar Sitamarhi Siwan Supaul Vaishali
Araria              0     0          0       0         0     0     -1        0
Arwal               0     0          0       0         0     0      0        0
Aurangabad          0     0          0       0         0     0      0        0
Banka               0     0          0       0         0     0      0        0
Begusarai          -1     0          0       0         0     0      0        0
Bhagalpur           0     0          0       0         0     0      0        0
           Balia Chandauli Deoria Gazipur Gorakhpur Kushinagar Maharajganj Mau
Araria         0         0      0       0         0          0           0   0
Arwal          0         0      0       0         0          0           0   0
Aurangabad     0         0      0       0         0          0           0   0
Banka          0         0      0       0         0          0           0   0
Begusarai      0         0      0       0         0          0           0   0
Bhagalpur      0         0      0       0         0          0           0   0
           Siddharthnagar
Araria                  0
Arwal                   0
Aurangabad              0
Banka                   0
Begusarai               0
Bhagalpur               0
Code 
## Also need to transform to factor for
## setting up the MRF smooth.
Irrig_Rev_rice_wheat$District <- as.factor(Irrig_Rev_rice_wheat$a_q103_district)

## Now note that not all regions are observed,
## therefore we need to remove those regions
## from the penalty matrix
rn <- rownames(K)
lv <- levels(Irrig_Rev_rice_wheat$District)
i <- rn %in% lv
K <- K[i, i]

set.seed(321)
library(bamlss)
b_rice_yield_MRF <- bamlss(f_rice_yield_MRF, data = Irrig_Rev_rice_wheat, family = "gaussian")
AICc 12433.92 logPost -6216.75 logLik -6063.08 edf 148.79 eps 2.3399 iteration   1
AICc 12239.92 logPost -5910.29 logLik -5967.69 edf 147.28 eps 0.4613 iteration   2
AICc 12211.69 logPost -5884.01 logLik -5956.93 edf 144.14 eps 0.1461 iteration   3
AICc 12205.11 logPost -5882.70 logLik -5953.35 edf 144.42 eps 2.9745 iteration   4
AICc 12203.12 logPost -5882.46 logLik -5952.28 edf 144.48 eps 0.0429 iteration   5
AICc 12202.50 logPost -5882.43 logLik -5951.95 edf 144.50 eps 0.0113 iteration   6
AICc 12202.32 logPost -5882.45 logLik -5951.84 edf 144.52 eps 0.0088 iteration   7
AICc 12202.28 logPost -5882.44 logLik -5951.82 edf 144.52 eps 0.0010 iteration   8
AICc 12202.27 logPost -5882.44 logLik -5951.81 edf 144.53 eps 0.0002 iteration   9
AICc 12202.27 logPost -5882.44 logLik -5951.81 edf 144.53 eps 0.0001 iteration  10
AICc 12202.27 logPost -5882.44 logLik -5951.81 edf 144.53 eps 0.0001 iteration  11
AICc 12202.26 logPost -5882.44 logLik -5951.80 edf 144.53 eps 0.0000 iteration  12
AICc 12202.26 logPost -5882.44 logLik -5951.80 edf 144.53 eps 0.0000 iteration  12
elapsed time: 10.33sec
Starting the sampler...

|                    |   0%  2.36min
|*                   |   5%  1.99min  6.28sec
|**                  |  10%  1.87min 12.47sec
|***                 |  15%  1.77min 18.75sec
|****                |  20%  1.67min 25.12sec
|*****               |  25%  1.58min 31.51sec
|******              |  30%  1.48min 38.17sec
|*******             |  35%  1.39min 44.97sec
|********            |  40%  1.30min 51.81sec
|*********           |  45%  1.19min 58.39sec
|**********          |  50%  1.08min  1.08min
|***********         |  55% 57.73sec  1.18min
|************        |  60% 51.27sec  1.28min
|*************       |  65% 45.03sec  1.39min
|**************      |  70% 38.91sec  1.51min
|***************     |  75% 32.49sec  1.62min
|****************    |  80% 25.98sec  1.73min
|*****************   |  85% 19.46sec  1.84min
|******************  |  90% 12.97sec  1.94min
|******************* |  95%  6.48sec  2.05min
|********************| 100%  0.00sec  2.16min
Code 
## First, note that we have the structured id = 'mrf1' and unstructured
## spatial effect id = 're2', also indicated in the model summary
summary(b_rice_yield_MRF)

Call:
bamlss(formula = f_rice_yield_MRF, family = "gaussian", data = Irrig_Rev_rice_wheat)
---
Family: gaussian 
Link function: mu = identity, sigma = log
*---
Formula mu:
---
b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + 
    s(sowdate_fmt_rice_day) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + 
    s(p2o5perha_rice) + s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
                            Mean    2.5%     50%   97.5% parameters
(Intercept)              1.90667 0.53173 1.64354 3.75448      0.013
rice_duration_class_long 0.07459 0.01275 0.07417 0.13858      0.075
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                       Mean      2.5%       50%     97.5%
s(sowdate_fmt_rice_day).tau21     5.545e-01 1.451e-02 3.038e-01 2.383e+00
s(sowdate_fmt_rice_day).alpha     1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(sowdate_fmt_rice_day).edf       3.793e+00 1.622e+00 3.716e+00 5.905e+00
s(g_q5305_irrig_times_rice).tau21 2.181e+00 2.788e-01 1.553e+00 7.637e+00
s(g_q5305_irrig_times_rice).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(g_q5305_irrig_times_rice).edf   6.445e+00 4.692e+00 6.475e+00 7.904e+00
s(nperha_rice).tau21              3.860e-01 1.403e-04 7.020e-02 2.928e+00
s(nperha_rice).alpha              1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(nperha_rice).edf                2.859e+00 1.009e+00 2.477e+00 6.347e+00
s(p2o5perha_rice).tau21           8.984e-02 8.439e-05 6.835e-03 8.142e-01
s(p2o5perha_rice).alpha           1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(p2o5perha_rice).edf             1.871e+00 1.007e+00 1.387e+00 4.623e+00
s(District,id='mrf1').tau21       3.920e-02 7.163e-05 1.029e-02 1.904e-01
s(District,id='mrf1').alpha       1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='mrf1').edf         2.229e+01 1.892e+00 2.823e+01 3.405e+01
s(District,id='re2').tau21        5.661e+00 2.042e-01 5.307e+00 1.505e+01
s(District,id='re2').alpha        1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='re2').edf          3.560e+01 3.390e+01 3.581e+01 3.594e+01
                                  parameters
s(sowdate_fmt_rice_day).tau21          0.182
s(sowdate_fmt_rice_day).alpha             NA
s(sowdate_fmt_rice_day).edf            3.291
s(g_q5305_irrig_times_rice).tau21      2.895
s(g_q5305_irrig_times_rice).alpha         NA
s(g_q5305_irrig_times_rice).edf        7.106
s(nperha_rice).tau21                   1.398
s(nperha_rice).alpha                      NA
s(nperha_rice).edf                     5.489
s(p2o5perha_rice).tau21                0.018
s(p2o5perha_rice).alpha                   NA
s(p2o5perha_rice).edf                  1.787
s(District,id='mrf1').tau21            0.029
s(District,id='mrf1').alpha               NA
s(District,id='mrf1').edf             31.761
s(District,id='re2').tau21            12.212
s(District,id='re2').alpha                NA
s(District,id='re2').edf              35.917
---
Formula sigma:
---
sigma ~ 1 + rice_duration_class_long + s(sowdate_fmt_rice_day) + 
    s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice) + 
    s(District, bs = "mrf", xt = list(penalty = K)) + s(District, 
    bs = "re")
-
Parametric coefficients:
                             Mean     2.5%      50%    97.5% parameters
(Intercept)              -0.12192 -0.17557 -0.12140 -0.06738     -0.114
rice_duration_class_long  0.06198  0.00957  0.06136  0.11587      0.060
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9858 0.8931 1.0000     1
-
Smooth terms:
                                       Mean      2.5%       50%     97.5%
s(sowdate_fmt_rice_day).tau21     2.045e-02 9.510e-05 3.988e-03 1.440e-01
s(sowdate_fmt_rice_day).alpha     9.779e-01 8.335e-01 9.997e-01 1.000e+00
s(sowdate_fmt_rice_day).edf       1.592e+00 1.009e+00 1.316e+00 3.418e+00
s(g_q5305_irrig_times_rice).tau21 2.123e-01 2.913e-04 1.141e-01 1.032e+00
s(g_q5305_irrig_times_rice).alpha 9.380e-01 6.284e-01 9.965e-01 1.000e+00
s(g_q5305_irrig_times_rice).edf   4.115e+00 1.138e+00 4.255e+00 6.535e+00
s(nperha_rice).tau21              1.830e-01 3.596e-04 3.058e-02 1.301e+00
s(nperha_rice).alpha              9.539e-01 6.658e-01 9.970e-01 1.000e+00
s(nperha_rice).edf                2.649e+00 1.037e+00 2.302e+00 6.011e+00
s(p2o5perha_rice).tau21           3.649e-01 1.628e-04 2.161e-02 2.782e+00
s(p2o5perha_rice).alpha           9.483e-01 5.766e-01 9.984e-01 1.000e+00
s(p2o5perha_rice).edf             2.717e+00 1.020e+00 2.109e+00 6.516e+00
s(District,id='mrf1').tau21       3.070e-03 8.543e-05 2.626e-03 8.938e-03
s(District,id='mrf1').alpha       7.837e-01 2.297e-01 8.795e-01 1.000e+00
s(District,id='mrf1').edf         2.155e+01 3.270e+00 2.360e+01 2.985e+01
s(District,id='re2').tau21        1.977e-02 3.995e-03 1.877e-02 4.442e-02
s(District,id='re2').alpha        7.394e-01 1.688e-01 8.059e-01 1.000e+00
s(District,id='re2').edf          2.747e+01 1.734e+01 2.857e+01 3.198e+01
                                  parameters
s(sowdate_fmt_rice_day).tau21          0.014
s(sowdate_fmt_rice_day).alpha             NA
s(sowdate_fmt_rice_day).edf            1.802
s(g_q5305_irrig_times_rice).tau21      0.191
s(g_q5305_irrig_times_rice).alpha         NA
s(g_q5305_irrig_times_rice).edf        4.768
s(nperha_rice).tau21                   0.091
s(nperha_rice).alpha                      NA
s(nperha_rice).edf                     3.139
s(p2o5perha_rice).tau21                2.253
s(p2o5perha_rice).alpha                   NA
s(p2o5perha_rice).edf                  6.288
s(District,id='mrf1').tau21            0.007
s(District,id='mrf1').alpha               NA
s(District,id='mrf1').edf             28.648
s(District,id='re2').tau21             0.002
s(District,id='re2').alpha                NA
s(District,id='re2').edf              10.540
---
Sampler summary:
-
DIC = 12123.09 logLik = -6011.782 pd = 99.5225
runtime = 130.58
---
Optimizer summary:
-
AICc = 12202.27 edf = 144.5353 logLik = -5951.809
logPost = -5882.444 nobs = 4537 runtime = 10.33
Code 
# Plot the nonlinear effect
plot(b_rice_yield_MRF, model = "mu", term = "s(sowdate_fmt_rice_day)")

Code 
plot(b_rice_yield_MRF, model = "mu", term = "s(g_q5305_irrig_times_rice)")

Code 
plot(b_rice_yield_MRF, model = "mu", term = "s(nperha_rice)")

Code 
plot(b_rice_yield_MRF, model = "mu", term = "s(p2o5perha_rice)")

2.2 Wheat model

Code 
f_wheat_yield_MRF <- list(
    l_ton_per_hectare ~ 1 + variety_type_NMWV + s(sowdate_fmt_wheat_day) + g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    sigma ~ 1 + variety_type_NMWV + s(sowdate_fmt_wheat_day) + g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re")
)
set.seed(321)
b_wheat_yield_MRF <- bamlss(f_wheat_yield_MRF, data = Irrig_Rev_rice_wheat, family = "gaussian")
AICc 8958.541 logPost -4419.70 logLik -4335.18 edf 139.61 eps 0.3246 iteration   1
AICc 8445.144 logPost -4096.84 logLik -4071.93 edf 145.76 eps 0.1431 iteration   2
AICc 8344.698 logPost -4042.08 logLik -4030.62 edf 137.40 eps 0.0345 iteration   3
AICc 8319.455 logPost -4032.32 logLik -4022.32 edf 133.33 eps 0.0132 iteration   4
AICc 8313.001 logPost -4030.95 logLik -4020.60 edf 131.91 eps 0.0050 iteration   5
AICc 8312.046 logPost -4029.90 logLik -4020.23 edf 131.81 eps 0.0019 iteration   6
AICc 8311.784 logPost -4029.54 logLik -4020.15 edf 131.76 eps 0.0010 iteration   7
AICc 8311.335 logPost -4029.64 logLik -4020.12 edf 131.57 eps 0.0005 iteration   8
AICc 8311.297 logPost -4029.53 logLik -4020.11 edf 131.57 eps 0.0003 iteration   9
AICc 8311.268 logPost -4029.43 logLik -4020.10 edf 131.56 eps 0.0002 iteration  10
AICc 8311.245 logPost -4029.35 logLik -4020.09 edf 131.56 eps 0.0002 iteration  11
AICc 8311.226 logPost -4029.28 logLik -4020.09 edf 131.56 eps 0.0002 iteration  12
AICc 8311.211 logPost -4029.22 logLik -4020.08 edf 131.56 eps 0.0001 iteration  13
AICc 8311.197 logPost -4029.16 logLik -4020.08 edf 131.55 eps 0.0001 iteration  14
AICc 8311.185 logPost -4029.11 logLik -4020.07 edf 131.55 eps 0.0001 iteration  15
AICc 8311.175 logPost -4029.06 logLik -4020.07 edf 131.55 eps 0.0001 iteration  16
AICc 8311.175 logPost -4029.06 logLik -4020.07 edf 131.55 eps 0.0001 iteration  16
elapsed time:  7.46sec
Starting the sampler...

|                    |   0%  1.41min
|*                   |   5%  1.38min  4.35sec
|**                  |  10%  1.28min  8.56sec
|***                 |  15%  1.24min 13.08sec
|****                |  20%  1.19min 17.89sec
|*****               |  25%  1.14min 22.74sec
|******              |  30%  1.07min 27.60sec
|*******             |  35%  1.00min 32.41sec
|********            |  40% 55.74sec 37.16sec
|*********           |  45% 51.26sec 41.94sec
|**********          |  50% 46.69sec 46.69sec
|***********         |  55% 41.91sec 51.22sec
|************        |  60% 37.43sec 56.14sec
|*************       |  65% 32.85sec  1.02min
|**************      |  70% 28.20sec  1.10min
|***************     |  75% 23.57sec  1.18min
|****************    |  80% 18.89sec  1.26min
|*****************   |  85% 14.14sec  1.34min
|******************  |  90%  9.43sec  1.41min
|******************* |  95%  4.71sec  1.49min
|********************| 100%  0.00sec  1.56min
Code 
## First, note that we have the structured id = 'mrf1' and unstructured
## spatial effect id = 're2', also indicated in the model summary
summary(b_wheat_yield_MRF)

Call:
bamlss(formula = f_wheat_yield_MRF, family = "gaussian", data = Irrig_Rev_rice_wheat)
---
Family: gaussian 
Link function: mu = identity, sigma = log
*---
Formula mu:
---
l_ton_per_hectare ~ 1 + variety_type_NMWV + s(sowdate_fmt_wheat_day) + 
    g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", 
    xt = list(penalty = K)) + s(District, bs = "re")
-
Parametric coefficients:
                          Mean       2.5%        50%      97.5% parameters
(Intercept)         -0.5717029 -1.1349600 -0.6585627  0.1644311     -1.359
variety_type_NMWV    0.3056401  0.2611485  0.3051494  0.3528470      0.303
g_q5305_irrig_times  0.4153518  0.3891090  0.4153159  0.4440488      0.411
nperha               0.0012344  0.0006667  0.0012304  0.0017908      0.001
p2o5perha            0.0031322  0.0020531  0.0031297  0.0042256      0.003
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                    Mean      2.5%       50%     97.5%
s(sowdate_fmt_wheat_day).tau21 4.212e-02 7.180e-05 3.403e-03 4.219e-01
s(sowdate_fmt_wheat_day).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(sowdate_fmt_wheat_day).edf   1.905e+00 1.014e+00 1.463e+00 5.031e+00
s(District,id='mrf1').tau21    3.313e-03 6.663e-05 9.203e-04 1.500e-02
s(District,id='mrf1').alpha    1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='mrf1').edf      1.869e+01 3.860e+00 1.867e+01 3.204e+01
s(District,id='re2').tau21     4.782e+00 1.753e+00 4.520e+00 9.634e+00
s(District,id='re2').alpha     1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='re2').edf       3.589e+01 3.576e+01 3.590e+01 3.596e+01
                               parameters
s(sowdate_fmt_wheat_day).tau21     15.401
s(sowdate_fmt_wheat_day).alpha         NA
s(sowdate_fmt_wheat_day).edf        8.309
s(District,id='mrf1').tau21         0.003
s(District,id='mrf1').alpha            NA
s(District,id='mrf1').edf          26.415
s(District,id='re2').tau21          6.085
s(District,id='re2').alpha             NA
s(District,id='re2').edf           35.897
---
Formula sigma:
---
sigma ~ 1 + variety_type_NMWV + s(sowdate_fmt_wheat_day) + g_q5305_irrig_times + 
    nperha + p2o5perha + s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
                          Mean       2.5%        50%      97.5% parameters
(Intercept)         -0.7813986 -0.8943190 -0.7812900 -0.6677184     -0.316
variety_type_NMWV    0.2055025  0.1500360  0.2056823  0.2547370      0.198
g_q5305_irrig_times  0.0924131  0.0622570  0.0926335  0.1238243      0.098
nperha              -0.0011585 -0.0018953 -0.0011631 -0.0004485     -0.001
p2o5perha            0.0017856  0.0003986  0.0017755  0.0031647      0.002
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.8959 0.4813 0.9768     1
-
Smooth terms:
                                    Mean      2.5%       50%     97.5%
s(sowdate_fmt_wheat_day).tau21 7.635e-02 6.028e-04 3.936e-02 3.405e-01
s(sowdate_fmt_wheat_day).alpha 9.637e-01 7.563e-01 9.985e-01 1.000e+00
s(sowdate_fmt_wheat_day).edf   2.514e+00 1.069e+00 2.487e+00 4.319e+00
s(District,id='mrf1').tau21    6.978e-04 5.086e-05 3.564e-04 2.818e-03
s(District,id='mrf1').alpha    8.951e-01 4.254e-01 9.764e-01 1.000e+00
s(District,id='mrf1').edf      1.105e+01 2.059e+00 9.448e+00 2.405e+01
s(District,id='re2').tau21     1.993e-02 9.130e-03 1.889e-02 3.532e-02
s(District,id='re2').alpha     7.365e-01 1.864e-01 8.248e-01 1.000e+00
s(District,id='re2').edf       2.836e+01 2.403e+01 2.860e+01 3.126e+01
                               parameters
s(sowdate_fmt_wheat_day).tau21      0.035
s(sowdate_fmt_wheat_day).alpha         NA
s(sowdate_fmt_wheat_day).edf        2.415
s(District,id='mrf1').tau21         0.001
s(District,id='mrf1').alpha            NA
s(District,id='mrf1').edf          13.967
s(District,id='re2').tau21          0.203
s(District,id='re2').alpha             NA
s(District,id='re2').edf           34.551
---
Sampler summary:
-
DIC = 8219.803 logLik = -4070.961 pd = 77.8806
runtime = 94.22
---
Optimizer summary:
-
AICc = 8311.175 edf = 131.5535 logLik = -4020.075
logPost = -4029.064 nobs = 4537 runtime = 7.46
Code 
# Plot the nonlinear effect
plot(b_wheat_yield_MRF, model = "mu", term = "s(sowdate_fmt_wheat_day)")

Code 
## Now, to predict the spatial effects we set up new data.
nd <- data.frame("District" = unique(Irrig_Rev_rice_wheat$District))

## Predict for the structured spatial effects.
p_str_wheat_yield_MRF <- predict(b_wheat_yield_MRF, newdata = nd, term = "s(District,id='mrf1')", intercept = FALSE)

## And the unstructured spatial effect.
p_unstr_wheat_yield_MRF <- predict(b_wheat_yield_MRF, newdata = nd, term = "s(District,id='re2')", intercept = FALSE)

## MRF smooth effect.
plotmap(India_aoi_sp_bnd,
    x = p_str_wheat_yield_MRF$mu, id = nd$District,
    main = expression(mu), title = "Structured spatial effect"
)

Code 
plotmap(India_aoi_sp_bnd,
    x = p_str_wheat_yield_MRF$sigma, id = nd$District,
    main = expression(sigma), title = "Structured spatial effect"
)

Code 
# Plotting using bivariate map
library(sf)
India_aoi_sf <- st_read("shp/India_aoi_sf_sp.shp")
Reading layer `India_aoi_sf_sp' from data source 
  `C:\Users\MMKONDIWA\OneDrive - CIMMYT\Documents\GitHub\Rice_Wheat_System_SDS_Tool\shp\India_aoi_sf_sp.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 47 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 82.42947 ymin: 24.28704 xmax: 88.29192 ymax: 27.85072
Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs
Code 
p_str_wheat_yield_MRF_dt <- as.data.frame(p_str_wheat_yield_MRF)

lv_dist <- as.data.frame(lv)
lv_dist_str <- cbind(lv_dist, p_str_wheat_yield_MRF_dt)


India_aoi_sf_selected_str <- merge(India_aoi_sf, lv_dist_str, by.x = "NAME_2", by.y = "lv")

library(biscale)
library(ggplot2)
data <- bi_class(India_aoi_sf_selected_str, x = "mu", y = "sigma", style = "quantile", dim = 3)
table(data$bi_class)

1-1 1-2 1-3 2-1 2-2 2-3 3-1 3-2 3-3 
  3   3   4   1   5   4   6   2   2
Code 
labels1 <- biscale::bi_class_breaks(
    data,
    x = mu,
    y = sigma,
    style = "quantile",
    dim = 3, dig_lab = 0, split = FALSE
)
labels1
$bi_x
[1] "-0.1--0.009" "-0.009-0.02" "0.02-0.1"   

$bi_y
[1] "-0.1--0.006"  "-0.006-0.008" "0.008-0.03"
Code 
previous_theme <- theme_set(theme_bw())

# create map
map <- ggplot() +
    geom_sf(data = data, mapping = aes(fill = bi_class), color = NA, size = 0.1, show.legend = FALSE) +
    bi_scale_fill(pal = "GrPink", dim = 3) +
    labs(title = "Structured spatial effect")
bi_theme()
List of 136
 $ line                            :List of 6
  ..$ colour       : chr "black"
  ..$ linewidth    : num 1.09
  ..$ linetype     : num 1
  ..$ lineend      : chr "butt"
  ..$ arrow        : logi FALSE
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_line" "element"
 $ rect                            :List of 5
  ..$ fill         : chr "white"
  ..$ colour       : chr "black"
  ..$ linewidth    : num 1.09
  ..$ linetype     : num 1
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_rect" "element"
 $ text                            :List of 11
  ..$ family       : chr "sans"
  ..$ face         : chr "plain"
  ..$ colour       : chr "#000000"
  ..$ size         : num 24
  ..$ hjust        : num 0.5
  ..$ vjust        : num 0.5
  ..$ angle        : num 0
  ..$ lineheight   : num 0.9
  ..$ margin       : 'margin' num [1:4] 0points 0points 0points 0points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : logi FALSE
  ..$ inherit.blank: logi FALSE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ title                           : NULL
 $ aspect.ratio                    : NULL
 $ axis.title                      : NULL
 $ axis.title.x                    :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : NULL
  ..$ vjust        : num 1
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 6points 0points 0points 0points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.title.x.top                :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : NULL
  ..$ vjust        : num 0
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 0points 0points 6points 0points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.title.x.bottom             : NULL
 $ axis.title.y                    :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : NULL
  ..$ vjust        : num 1
  ..$ angle        : num 90
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 0points 6points 0points 0points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.title.y.left               : NULL
 $ axis.title.y.right              :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : NULL
  ..$ vjust        : num 1
  ..$ angle        : num -90
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 0points 0points 0points 6points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.text                       : list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 $ axis.text.x                     :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : NULL
  ..$ vjust        : num 1
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 4.8points 0points 0points 0points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.text.x.top                 :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : NULL
  ..$ vjust        : num 0
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 0points 0points 4.8points 0points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.text.x.bottom              : NULL
 $ axis.text.y                     :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : num 1
  ..$ vjust        : NULL
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 0points 4.8points 0points 0points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.text.y.left                : NULL
 $ axis.text.y.right               :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : num 0
  ..$ vjust        : NULL
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 0points 0points 0points 4.8points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.text.theta                 : NULL
 $ axis.text.r                     :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : num 0.5
  ..$ vjust        : NULL
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 0points 4.8points 0points 4.8points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.ticks                      : list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 $ axis.ticks.x                    : NULL
 $ axis.ticks.x.top                : NULL
 $ axis.ticks.x.bottom             : NULL
 $ axis.ticks.y                    : NULL
 $ axis.ticks.y.left               : NULL
 $ axis.ticks.y.right              : NULL
 $ axis.ticks.theta                : NULL
 $ axis.ticks.r                    : NULL
 $ axis.minor.ticks.x.top          : NULL
 $ axis.minor.ticks.x.bottom       : NULL
 $ axis.minor.ticks.y.left         : NULL
 $ axis.minor.ticks.y.right        : NULL
 $ axis.minor.ticks.theta          : NULL
 $ axis.minor.ticks.r              : NULL
 $ axis.ticks.length               : 'simpleUnit' num 6points
  ..- attr(*, "unit")= int 8
 $ axis.ticks.length.x             : NULL
 $ axis.ticks.length.x.top         : NULL
 $ axis.ticks.length.x.bottom      : NULL
 $ axis.ticks.length.y             : NULL
 $ axis.ticks.length.y.left        : NULL
 $ axis.ticks.length.y.right       : NULL
 $ axis.ticks.length.theta         : NULL
 $ axis.ticks.length.r             : NULL
 $ axis.minor.ticks.length         : 'rel' num 0.75
 $ axis.minor.ticks.length.x       : NULL
 $ axis.minor.ticks.length.x.top   : NULL
 $ axis.minor.ticks.length.x.bottom: NULL
 $ axis.minor.ticks.length.y       : NULL
 $ axis.minor.ticks.length.y.left  : NULL
 $ axis.minor.ticks.length.y.right : NULL
 $ axis.minor.ticks.length.theta   : NULL
 $ axis.minor.ticks.length.r       : NULL
 $ axis.line                       : list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 $ axis.line.x                     : NULL
 $ axis.line.x.top                 : NULL
 $ axis.line.x.bottom              : NULL
 $ axis.line.y                     : NULL
 $ axis.line.y.left                : NULL
 $ axis.line.y.right               : NULL
 $ axis.line.theta                 : NULL
 $ axis.line.r                     : NULL
 $ legend.background               :List of 5
  ..$ fill         : chr "#ffffff"
  ..$ colour       : logi NA
  ..$ linewidth    : NULL
  ..$ linetype     : NULL
  ..$ inherit.blank: logi FALSE
  ..- attr(*, "class")= chr [1:2] "element_rect" "element"
 $ legend.margin                   : 'margin' num [1:4] 12points 12points 12points 12points
  ..- attr(*, "unit")= int 8
 $ legend.spacing                  : 'simpleUnit' num 24points
  ..- attr(*, "unit")= int 8
 $ legend.spacing.x                : NULL
 $ legend.spacing.y                : NULL
 $ legend.key                      : list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 $ legend.key.size                 : 'simpleUnit' num 1.2lines
  ..- attr(*, "unit")= int 3
 $ legend.key.height               : NULL
 $ legend.key.width                : NULL
 $ legend.key.spacing              : 'simpleUnit' num 12points
  ..- attr(*, "unit")= int 8
 $ legend.key.spacing.x            : NULL
 $ legend.key.spacing.y            : NULL
 $ legend.frame                    : NULL
 $ legend.ticks                    : NULL
 $ legend.ticks.length             : 'rel' num 0.2
 $ legend.axis.line                : NULL
 $ legend.text                     :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : chr "#000000"
  ..$ size         : 'rel' num 0.8
  ..$ hjust        : num 0
  ..$ vjust        : NULL
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : NULL
  ..$ debug        : NULL
  ..$ inherit.blank: logi FALSE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ legend.text.position            : NULL
 $ legend.title                    :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : chr "#000000"
  ..$ size         : NULL
  ..$ hjust        : num 0
  ..$ vjust        : NULL
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : NULL
  ..$ debug        : NULL
  ..$ inherit.blank: logi FALSE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ legend.title.position           : NULL
 $ legend.position                 : chr "right"
 $ legend.position.inside          : NULL
 $ legend.direction                : NULL
 $ legend.byrow                    : NULL
 $ legend.justification            : chr "center"
 $ legend.justification.top        : NULL
 $ legend.justification.bottom     : NULL
 $ legend.justification.left       : NULL
 $ legend.justification.right      : NULL
 $ legend.justification.inside     : NULL
 $ legend.location                 : NULL
 $ legend.box                      : NULL
 $ legend.box.just                 : NULL
 $ legend.box.margin               : 'margin' num [1:4] 0cm 0cm 0cm 0cm
  ..- attr(*, "unit")= int 1
 $ legend.box.background           : list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 $ legend.box.spacing              : 'simpleUnit' num 24points
  ..- attr(*, "unit")= int 8
  [list output truncated]
 - attr(*, "class")= chr [1:2] "theme" "gg"
 - attr(*, "complete")= logi TRUE
 - attr(*, "validate")= logi TRUE
Code 
map

Code 
legend <- bi_legend(
    pal = "GrPink",
    breaks = labels1,
    xlab = "Higher mu",
    ylab = "Higher sigma",
    size = 12
)
legend

Code 
# combine map with legend
library(cowplot)
finalPlot <- ggdraw() +
    draw_plot(map, 0, 0, 1, 1) +
    draw_plot(legend, 0.7, 0.7, 0.2, 0.2)

finalPlot

Code 
## Random effects.
plotmap(India_aoi_sp_bnd,
    x = p_unstr_wheat_yield_MRF$mu, id = nd$District, title = "Unstructured spatial effect"
)

Code 
plotmap(India_aoi_sp_bnd,
    x = p_unstr_wheat_yield_MRF$sigma, id = nd$District, title = "Unstructured spatial effect"
)

Code 
# Plotting using bivariate map
library(sf)
India_aoi_sf <- st_read("shp/India_aoi_sf_sp.shp")
Reading layer `India_aoi_sf_sp' from data source 
  `C:\Users\MMKONDIWA\OneDrive - CIMMYT\Documents\GitHub\Rice_Wheat_System_SDS_Tool\shp\India_aoi_sf_sp.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 47 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 82.42947 ymin: 24.28704 xmax: 88.29192 ymax: 27.85072
Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs
Code 
p_unstr_wheat_yield_MRF_dt <- as.data.frame(p_unstr_wheat_yield_MRF)

lv_dist <- as.data.frame(lv)
lv_dist_unstr <- cbind(lv_dist, p_unstr_wheat_yield_MRF_dt)

library(dplyr)
lv_dist_unstr <- rename(lv_dist_unstr, mu_unstr = mu)
lv_dist_unstr <- rename(lv_dist_unstr, sigma_unstr = sigma)

India_aoi_sf_selected_unstr <- merge(India_aoi_sf, lv_dist_unstr, by.x = "NAME_2", by.y = "lv")

library(biscale)
library(ggplot2)
data <- bi_class(India_aoi_sf_selected_unstr, x = "mu_unstr", y = "sigma_unstr", style = "quantile", dim = 3)
table(data$bi_class)

1-1 1-2 1-3 2-1 2-2 2-3 3-1 3-2 3-3 
  4   4   2   3   3   4   3   3   4
Code 
previous_theme <- theme_set(theme_bw())

# create map
map <- ggplot() +
    geom_sf(data = data, mapping = aes(fill = bi_class), color = NA, size = 0.1, show.legend = FALSE) +
    bi_scale_fill(pal = "GrPink", dim = 3) +
    labs(title = "Unstructured spatial effect")
bi_theme()
List of 136
 $ line                            :List of 6
  ..$ colour       : chr "black"
  ..$ linewidth    : num 1.09
  ..$ linetype     : num 1
  ..$ lineend      : chr "butt"
  ..$ arrow        : logi FALSE
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_line" "element"
 $ rect                            :List of 5
  ..$ fill         : chr "white"
  ..$ colour       : chr "black"
  ..$ linewidth    : num 1.09
  ..$ linetype     : num 1
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_rect" "element"
 $ text                            :List of 11
  ..$ family       : chr "sans"
  ..$ face         : chr "plain"
  ..$ colour       : chr "#000000"
  ..$ size         : num 24
  ..$ hjust        : num 0.5
  ..$ vjust        : num 0.5
  ..$ angle        : num 0
  ..$ lineheight   : num 0.9
  ..$ margin       : 'margin' num [1:4] 0points 0points 0points 0points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : logi FALSE
  ..$ inherit.blank: logi FALSE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ title                           : NULL
 $ aspect.ratio                    : NULL
 $ axis.title                      : NULL
 $ axis.title.x                    :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : NULL
  ..$ vjust        : num 1
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 6points 0points 0points 0points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.title.x.top                :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : NULL
  ..$ vjust        : num 0
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 0points 0points 6points 0points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.title.x.bottom             : NULL
 $ axis.title.y                    :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : NULL
  ..$ vjust        : num 1
  ..$ angle        : num 90
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 0points 6points 0points 0points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.title.y.left               : NULL
 $ axis.title.y.right              :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : NULL
  ..$ vjust        : num 1
  ..$ angle        : num -90
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 0points 0points 0points 6points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.text                       : list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 $ axis.text.x                     :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : NULL
  ..$ vjust        : num 1
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 4.8points 0points 0points 0points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.text.x.top                 :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : NULL
  ..$ vjust        : num 0
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 0points 0points 4.8points 0points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.text.x.bottom              : NULL
 $ axis.text.y                     :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : num 1
  ..$ vjust        : NULL
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 0points 4.8points 0points 0points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.text.y.left                : NULL
 $ axis.text.y.right               :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : num 0
  ..$ vjust        : NULL
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 0points 0points 0points 4.8points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.text.theta                 : NULL
 $ axis.text.r                     :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : num 0.5
  ..$ vjust        : NULL
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 0points 4.8points 0points 4.8points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.ticks                      : list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 $ axis.ticks.x                    : NULL
 $ axis.ticks.x.top                : NULL
 $ axis.ticks.x.bottom             : NULL
 $ axis.ticks.y                    : NULL
 $ axis.ticks.y.left               : NULL
 $ axis.ticks.y.right              : NULL
 $ axis.ticks.theta                : NULL
 $ axis.ticks.r                    : NULL
 $ axis.minor.ticks.x.top          : NULL
 $ axis.minor.ticks.x.bottom       : NULL
 $ axis.minor.ticks.y.left         : NULL
 $ axis.minor.ticks.y.right        : NULL
 $ axis.minor.ticks.theta          : NULL
 $ axis.minor.ticks.r              : NULL
 $ axis.ticks.length               : 'simpleUnit' num 6points
  ..- attr(*, "unit")= int 8
 $ axis.ticks.length.x             : NULL
 $ axis.ticks.length.x.top         : NULL
 $ axis.ticks.length.x.bottom      : NULL
 $ axis.ticks.length.y             : NULL
 $ axis.ticks.length.y.left        : NULL
 $ axis.ticks.length.y.right       : NULL
 $ axis.ticks.length.theta         : NULL
 $ axis.ticks.length.r             : NULL
 $ axis.minor.ticks.length         : 'rel' num 0.75
 $ axis.minor.ticks.length.x       : NULL
 $ axis.minor.ticks.length.x.top   : NULL
 $ axis.minor.ticks.length.x.bottom: NULL
 $ axis.minor.ticks.length.y       : NULL
 $ axis.minor.ticks.length.y.left  : NULL
 $ axis.minor.ticks.length.y.right : NULL
 $ axis.minor.ticks.length.theta   : NULL
 $ axis.minor.ticks.length.r       : NULL
 $ axis.line                       : list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 $ axis.line.x                     : NULL
 $ axis.line.x.top                 : NULL
 $ axis.line.x.bottom              : NULL
 $ axis.line.y                     : NULL
 $ axis.line.y.left                : NULL
 $ axis.line.y.right               : NULL
 $ axis.line.theta                 : NULL
 $ axis.line.r                     : NULL
 $ legend.background               :List of 5
  ..$ fill         : chr "#ffffff"
  ..$ colour       : logi NA
  ..$ linewidth    : NULL
  ..$ linetype     : NULL
  ..$ inherit.blank: logi FALSE
  ..- attr(*, "class")= chr [1:2] "element_rect" "element"
 $ legend.margin                   : 'margin' num [1:4] 12points 12points 12points 12points
  ..- attr(*, "unit")= int 8
 $ legend.spacing                  : 'simpleUnit' num 24points
  ..- attr(*, "unit")= int 8
 $ legend.spacing.x                : NULL
 $ legend.spacing.y                : NULL
 $ legend.key                      : list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 $ legend.key.size                 : 'simpleUnit' num 1.2lines
  ..- attr(*, "unit")= int 3
 $ legend.key.height               : NULL
 $ legend.key.width                : NULL
 $ legend.key.spacing              : 'simpleUnit' num 12points
  ..- attr(*, "unit")= int 8
 $ legend.key.spacing.x            : NULL
 $ legend.key.spacing.y            : NULL
 $ legend.frame                    : NULL
 $ legend.ticks                    : NULL
 $ legend.ticks.length             : 'rel' num 0.2
 $ legend.axis.line                : NULL
 $ legend.text                     :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : chr "#000000"
  ..$ size         : 'rel' num 0.8
  ..$ hjust        : num 0
  ..$ vjust        : NULL
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : NULL
  ..$ debug        : NULL
  ..$ inherit.blank: logi FALSE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ legend.text.position            : NULL
 $ legend.title                    :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : chr "#000000"
  ..$ size         : NULL
  ..$ hjust        : num 0
  ..$ vjust        : NULL
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : NULL
  ..$ debug        : NULL
  ..$ inherit.blank: logi FALSE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ legend.title.position           : NULL
 $ legend.position                 : chr "right"
 $ legend.position.inside          : NULL
 $ legend.direction                : NULL
 $ legend.byrow                    : NULL
 $ legend.justification            : chr "center"
 $ legend.justification.top        : NULL
 $ legend.justification.bottom     : NULL
 $ legend.justification.left       : NULL
 $ legend.justification.right      : NULL
 $ legend.justification.inside     : NULL
 $ legend.location                 : NULL
 $ legend.box                      : NULL
 $ legend.box.just                 : NULL
 $ legend.box.margin               : 'margin' num [1:4] 0cm 0cm 0cm 0cm
  ..- attr(*, "unit")= int 1
 $ legend.box.background           : list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 $ legend.box.spacing              : 'simpleUnit' num 24points
  ..- attr(*, "unit")= int 8
  [list output truncated]
 - attr(*, "class")= chr [1:2] "theme" "gg"
 - attr(*, "complete")= logi TRUE
 - attr(*, "validate")= logi TRUE
Code 
map

Code 
labels1 <- biscale::bi_class_breaks(
    data,
    x = mu_unstr,
    y = sigma_unstr,
    style = "quantile",
    dim = 3, dig_lab = 0, split = FALSE
)
labels1
$bi_x
[1] "1.5-2"   "2-2.2"   "2.2-2.6"

$bi_y
[1] "-0.3--0.04" "-0.04-0.05" "0.05-0.1"
Code 
legend <- bi_legend(
    pal = "GrPink",
    breaks = labels1,
    dim = 3,
    xlab = "Higher mu",
    ylab = "Higher sigma",
    size = 10
)

# combine map with legend
library(cowplot)
finalPlot <- ggdraw() +
    draw_plot(map, 0, 0, 1, 1) +
    draw_plot(legend, 0.7, 0.7, 0.2, 0.2)

finalPlot

Code 
#

# Structured vs unstructured

India_aoi_sf_selected_str_dt <- subset(India_aoi_sf_selected_str, select = c("NAME_2", "mu", "sigma"))

India_aoi_sf_selected_str_dt$geometry <- NULL

India_aoi_sf_selected_unstr_str <- merge(India_aoi_sf_selected_unstr, India_aoi_sf_selected_str_dt, by = "NAME_2")

library(biscale)
library(ggplot2)
data <- bi_class(India_aoi_sf_selected_unstr_str, x = "mu", y = "mu_unstr", style = "quantile", dim = 3)
table(data$bi_class)

1-1 1-2 2-1 2-2 2-3 3-1 3-2 3-3 
  4   6   4   2   4   2   2   6
Code 
previous_theme <- theme_set(theme_bw())

# create map
map <- ggplot() +
    geom_sf(data = data, mapping = aes(fill = bi_class), color = NA, size = 0.1, show.legend = FALSE) +
    bi_scale_fill(pal = "GrPink", dim = 3) +
    labs(title = "Structured and unstructured spatial effect")
bi_theme()
List of 136
 $ line                            :List of 6
  ..$ colour       : chr "black"
  ..$ linewidth    : num 1.09
  ..$ linetype     : num 1
  ..$ lineend      : chr "butt"
  ..$ arrow        : logi FALSE
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_line" "element"
 $ rect                            :List of 5
  ..$ fill         : chr "white"
  ..$ colour       : chr "black"
  ..$ linewidth    : num 1.09
  ..$ linetype     : num 1
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_rect" "element"
 $ text                            :List of 11
  ..$ family       : chr "sans"
  ..$ face         : chr "plain"
  ..$ colour       : chr "#000000"
  ..$ size         : num 24
  ..$ hjust        : num 0.5
  ..$ vjust        : num 0.5
  ..$ angle        : num 0
  ..$ lineheight   : num 0.9
  ..$ margin       : 'margin' num [1:4] 0points 0points 0points 0points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : logi FALSE
  ..$ inherit.blank: logi FALSE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ title                           : NULL
 $ aspect.ratio                    : NULL
 $ axis.title                      : NULL
 $ axis.title.x                    :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : NULL
  ..$ vjust        : num 1
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 6points 0points 0points 0points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.title.x.top                :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : NULL
  ..$ vjust        : num 0
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 0points 0points 6points 0points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.title.x.bottom             : NULL
 $ axis.title.y                    :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : NULL
  ..$ vjust        : num 1
  ..$ angle        : num 90
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 0points 6points 0points 0points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.title.y.left               : NULL
 $ axis.title.y.right              :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : NULL
  ..$ vjust        : num 1
  ..$ angle        : num -90
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 0points 0points 0points 6points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.text                       : list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 $ axis.text.x                     :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : NULL
  ..$ vjust        : num 1
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 4.8points 0points 0points 0points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.text.x.top                 :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : NULL
  ..$ vjust        : num 0
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 0points 0points 4.8points 0points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.text.x.bottom              : NULL
 $ axis.text.y                     :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : num 1
  ..$ vjust        : NULL
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 0points 4.8points 0points 0points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.text.y.left                : NULL
 $ axis.text.y.right               :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : num 0
  ..$ vjust        : NULL
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 0points 0points 0points 4.8points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.text.theta                 : NULL
 $ axis.text.r                     :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : NULL
  ..$ size         : NULL
  ..$ hjust        : num 0.5
  ..$ vjust        : NULL
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : 'margin' num [1:4] 0points 4.8points 0points 4.8points
  .. ..- attr(*, "unit")= int 8
  ..$ debug        : NULL
  ..$ inherit.blank: logi TRUE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ axis.ticks                      : list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 $ axis.ticks.x                    : NULL
 $ axis.ticks.x.top                : NULL
 $ axis.ticks.x.bottom             : NULL
 $ axis.ticks.y                    : NULL
 $ axis.ticks.y.left               : NULL
 $ axis.ticks.y.right              : NULL
 $ axis.ticks.theta                : NULL
 $ axis.ticks.r                    : NULL
 $ axis.minor.ticks.x.top          : NULL
 $ axis.minor.ticks.x.bottom       : NULL
 $ axis.minor.ticks.y.left         : NULL
 $ axis.minor.ticks.y.right        : NULL
 $ axis.minor.ticks.theta          : NULL
 $ axis.minor.ticks.r              : NULL
 $ axis.ticks.length               : 'simpleUnit' num 6points
  ..- attr(*, "unit")= int 8
 $ axis.ticks.length.x             : NULL
 $ axis.ticks.length.x.top         : NULL
 $ axis.ticks.length.x.bottom      : NULL
 $ axis.ticks.length.y             : NULL
 $ axis.ticks.length.y.left        : NULL
 $ axis.ticks.length.y.right       : NULL
 $ axis.ticks.length.theta         : NULL
 $ axis.ticks.length.r             : NULL
 $ axis.minor.ticks.length         : 'rel' num 0.75
 $ axis.minor.ticks.length.x       : NULL
 $ axis.minor.ticks.length.x.top   : NULL
 $ axis.minor.ticks.length.x.bottom: NULL
 $ axis.minor.ticks.length.y       : NULL
 $ axis.minor.ticks.length.y.left  : NULL
 $ axis.minor.ticks.length.y.right : NULL
 $ axis.minor.ticks.length.theta   : NULL
 $ axis.minor.ticks.length.r       : NULL
 $ axis.line                       : list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 $ axis.line.x                     : NULL
 $ axis.line.x.top                 : NULL
 $ axis.line.x.bottom              : NULL
 $ axis.line.y                     : NULL
 $ axis.line.y.left                : NULL
 $ axis.line.y.right               : NULL
 $ axis.line.theta                 : NULL
 $ axis.line.r                     : NULL
 $ legend.background               :List of 5
  ..$ fill         : chr "#ffffff"
  ..$ colour       : logi NA
  ..$ linewidth    : NULL
  ..$ linetype     : NULL
  ..$ inherit.blank: logi FALSE
  ..- attr(*, "class")= chr [1:2] "element_rect" "element"
 $ legend.margin                   : 'margin' num [1:4] 12points 12points 12points 12points
  ..- attr(*, "unit")= int 8
 $ legend.spacing                  : 'simpleUnit' num 24points
  ..- attr(*, "unit")= int 8
 $ legend.spacing.x                : NULL
 $ legend.spacing.y                : NULL
 $ legend.key                      : list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 $ legend.key.size                 : 'simpleUnit' num 1.2lines
  ..- attr(*, "unit")= int 3
 $ legend.key.height               : NULL
 $ legend.key.width                : NULL
 $ legend.key.spacing              : 'simpleUnit' num 12points
  ..- attr(*, "unit")= int 8
 $ legend.key.spacing.x            : NULL
 $ legend.key.spacing.y            : NULL
 $ legend.frame                    : NULL
 $ legend.ticks                    : NULL
 $ legend.ticks.length             : 'rel' num 0.2
 $ legend.axis.line                : NULL
 $ legend.text                     :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : chr "#000000"
  ..$ size         : 'rel' num 0.8
  ..$ hjust        : num 0
  ..$ vjust        : NULL
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : NULL
  ..$ debug        : NULL
  ..$ inherit.blank: logi FALSE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ legend.text.position            : NULL
 $ legend.title                    :List of 11
  ..$ family       : NULL
  ..$ face         : NULL
  ..$ colour       : chr "#000000"
  ..$ size         : NULL
  ..$ hjust        : num 0
  ..$ vjust        : NULL
  ..$ angle        : NULL
  ..$ lineheight   : NULL
  ..$ margin       : NULL
  ..$ debug        : NULL
  ..$ inherit.blank: logi FALSE
  ..- attr(*, "class")= chr [1:2] "element_text" "element"
 $ legend.title.position           : NULL
 $ legend.position                 : chr "right"
 $ legend.position.inside          : NULL
 $ legend.direction                : NULL
 $ legend.byrow                    : NULL
 $ legend.justification            : chr "center"
 $ legend.justification.top        : NULL
 $ legend.justification.bottom     : NULL
 $ legend.justification.left       : NULL
 $ legend.justification.right      : NULL
 $ legend.justification.inside     : NULL
 $ legend.location                 : NULL
 $ legend.box                      : NULL
 $ legend.box.just                 : NULL
 $ legend.box.margin               : 'margin' num [1:4] 0cm 0cm 0cm 0cm
  ..- attr(*, "unit")= int 1
 $ legend.box.background           : list()
  ..- attr(*, "class")= chr [1:2] "element_blank" "element"
 $ legend.box.spacing              : 'simpleUnit' num 24points
  ..- attr(*, "unit")= int 8
  [list output truncated]
 - attr(*, "class")= chr [1:2] "theme" "gg"
 - attr(*, "complete")= logi TRUE
 - attr(*, "validate")= logi TRUE
Code 
map

Code 
labels1 <- biscale::bi_class_breaks(
    data,
    x = mu,
    y = mu_unstr,
    style = "quantile",
    dim = 3, dig_lab = 0, split = FALSE
)
labels1
$bi_x
[1] "-0.1--0.009" "-0.009-0.02" "0.02-0.1"   

$bi_y
[1] "1.5-2"   "2-2.2"   "2.2-2.6"
Code 
legend <- bi_legend(
    pal = "GrPink",
    breaks = labels1,
    dim = 3,
    xlab = "Higher mu str",
    ylab = "Higher mu unstr",
    size = 12
)

# combine map with legend
library(cowplot)
finalPlot <- ggdraw() +
    draw_plot(map, 0, 0, 1, 1) +
    draw_plot(legend, 0.7, 0.7, 0.2, 0.2)

finalPlot

3 Multivariate Spatial Geoadditive Bayesian Regression Model

4 Spatial + for multivariate model

Code 
library(sp)
library(mgcv)
library(bamlss)

# Multivariate geoadditive model
# remotes::install_git("https://git.uibk.ac.at/c4031039/mvnchol")
# install_github("https://github.com/meteosimon/mvnchol")
library(mvnchol)

# Remove NAs in the monsoon variables
Irrig_Rev_rice_wheat <- subset(Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$onset_2017)))
Irrig_Rev_rice_wheat <- subset(Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$monsoon_onset_dev)))
Irrig_Rev_rice_wheat <- subset(Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$median_onset_82_15)))
Irrig_Rev_rice_wheat <- subset(Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$sd_onset_82_15)))

# Rice first stage
f_sow_rice_1st_stage <- list(
    sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) + s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + s(sd_onset_82_15) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re")
)

f_sow_rice_1st_stage_MRF <- bamlss(f_sow_rice_1st_stage, data = Irrig_Rev_rice_wheat, family = "gaussian")
AICc 38862.03 logPost -165306. logLik -19371.2 edf 59.010 eps 0.1226 iteration   1
AICc 34841.08 logPost -31072.7 logLik -17339.3 edf 79.788 eps 0.0980 iteration   2
AICc 33498.94 logPost -18399.0 logLik -16648.2 edf 98.970 eps 0.0632 iteration   3
AICc 33180.98 logPost -17093.4 logLik -16491.0 edf 97.287 eps 0.0297 iteration   4
AICc 33096.83 logPost -16986.5 logLik -16452.6 edf 93.751 eps 0.0069 iteration   5
AICc 33056.53 logPost -16935.2 logLik -16435.7 edf 90.674 eps 0.0017 iteration   6
AICc 33040.37 logPost -16899.0 logLik -16429.8 edf 88.587 eps 0.0009 iteration   7
AICc 33035.63 logPost -16869.0 logLik -16428.1 edf 87.868 eps 0.0003 iteration   8
AICc 33028.91 logPost -16845.9 logLik -16427.6 edf 85.110 eps 0.0001 iteration   9
AICc 33028.60 logPost -16840.2 logLik -16427.5 edf 85.077 eps 0.0000 iteration  10
AICc 33028.60 logPost -16840.2 logLik -16427.5 edf 85.077 eps 0.0000 iteration  10
elapsed time:  4.35sec
Starting the sampler...

|                    |   0% 59.50sec
|*                   |   5% 59.66sec  3.14sec
|**                  |  10% 56.43sec  6.27sec
|***                 |  15% 52.36sec  9.24sec
|****                |  20% 51.32sec 12.83sec
|*****               |  25% 49.35sec 16.45sec
|******              |  30% 46.22sec 19.81sec
|*******             |  35% 43.29sec 23.31sec
|********            |  40% 40.11sec 26.74sec
|*********           |  45% 36.91sec 30.20sec
|**********          |  50% 33.63sec 33.63sec
|***********         |  55% 30.50sec 37.28sec
|************        |  60% 27.20sec 40.80sec
|*************       |  65% 23.90sec 44.39sec
|**************      |  70% 20.55sec 47.95sec
|***************     |  75% 17.19sec 51.58sec
|****************    |  80% 13.77sec 55.06sec
|*****************   |  85% 10.28sec 58.24sec
|******************  |  90%  6.84sec  1.03min
|******************* |  95%  3.40sec  1.08min
|********************| 100%  0.00sec  1.13min
Code 
fitted_f_sow_rice_1st_stage_MRF <- f_sow_rice_1st_stage_MRF$fitted

Irrig_Rev_rice_wheat$Res_rice_sow <- Irrig_Rev_rice_wheat$sowdate_fmt_rice_day - fitted_f_sow_rice_1st_stage_MRF$mu

# Wheat first stage
f_sow_wheat_1st_stage <- list(
    sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) + s(gw_2018) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re")
)

f_sow_wheat_1st_stage_MRF <- bamlss(f_sow_wheat_1st_stage, data = Irrig_Rev_rice_wheat, family = "gaussian")
AICc 49188.65 logPost -488850. logLik -24543.4 edf 50.285 eps 0.2566 iteration   1
AICc 44824.42 logPost -64385.9 logLik -22386.0 edf 26.046 eps 0.0753 iteration   2
AICc 40409.05 logPost -24539.6 logLik -20136.0 edf 67.444 eps 0.0824 iteration   3
AICc 36648.46 logPost -18973.1 logLik -18246.6 edf 76.219 eps 0.0789 iteration   4
AICc 33900.77 logPost -17313.2 logLik -16871.9 edf 77.104 eps 0.0796 iteration   5
AICc 32636.01 logPost -16736.4 logLik -16238.1 edf 78.490 eps 0.0660 iteration   6
AICc 32419.72 logPost -16629.9 logLik -16129.2 edf 79.163 eps 0.0324 iteration   7
AICc 32413.31 logPost -16627.4 logLik -16125.7 edf 79.502 eps 0.0056 iteration   8
AICc 32412.92 logPost -16627.3 logLik -16125.4 edf 79.574 eps 0.0001 iteration   9
AICc 32412.79 logPost -16627.2 logLik -16125.3 edf 79.581 eps 0.0000 iteration  10
AICc 32412.79 logPost -16627.2 logLik -16125.3 edf 79.581 eps 0.0000 iteration  10
elapsed time:  2.92sec
Starting the sampler...

|                    |   0% 44.03sec
|*                   |   5% 44.84sec  2.36sec
|**                  |  10% 43.38sec  4.82sec
|***                 |  15% 41.08sec  7.25sec
|****                |  20% 39.80sec  9.95sec
|*****               |  25% 38.25sec 12.75sec
|******              |  30% 36.52sec 15.65sec
|*******             |  35% 34.30sec 18.47sec
|********            |  40% 31.98sec 21.32sec
|*********           |  45% 29.21sec 23.90sec
|**********          |  50% 26.61sec 26.61sec
|***********         |  55% 24.16sec 29.53sec
|************        |  60% 21.78sec 32.67sec
|*************       |  65% 19.15sec 35.56sec
|**************      |  70% 16.49sec 38.48sec
|***************     |  75% 13.84sec 41.53sec
|****************    |  80% 11.07sec 44.28sec
|*****************   |  85%  8.28sec 46.90sec
|******************  |  90%  5.52sec 49.70sec
|******************* |  95%  2.75sec 52.32sec
|********************| 100%  0.00sec 54.93sec
Code 
fitted_f_sow_wheat_1st_stage_MRF <- f_sow_wheat_1st_stage_MRF$fitted

Irrig_Rev_rice_wheat$Res_wheat_sow <- Irrig_Rev_rice_wheat$sowdate_fmt_wheat_day - fitted_f_sow_wheat_1st_stage_MRF$mu


# Multivariate with sowing dates as endogenous
f_rice_wheat_yield_MRF_corr_endo <- list(
    sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) + s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + s(sd_onset_82_15) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    l_ton_per_hectare ~ 1 + variety_type_NMWV + s(Res_wheat_sow) + g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", xt = list("penalty" = K)) + s(District, bs = "re"),
    lamdiag1 ~ 1,
    lamdiag2 ~ 1,
    lamdiag3 ~ 1,
    lamdiag4 ~ 1,
    lambda12 ~ 1,
    lambda13 ~ 1,
    lambda14 ~ 1,
    lambda23 ~ 1,
    lambda24 ~ s(District, bs = "mrf", xt = list("penalty" = K)) + s(District, bs = "re"),
    lambda34 ~ 1
)

multivariate_geo_sow_MRF_corr_endo <- bamlss(f_rice_wheat_yield_MRF_corr_endo, type = "modified", family = mvnchol_bamlss(k = 4), data = Irrig_Rev_rice_wheat)
AICc 349476.1 logPost -2543606 logLik -174262. edf 430.14 eps 1.0000 iteration   1
AICc 157743.5 logPost -2447306 logLik -78382.3 edf 441.55 eps 0.1825 iteration   2
AICc 104307.5 logPost -2420173 logLik -51667.3 edf 439.10 eps 0.2070 iteration   3
AICc 89316.74 logPost -2412636 logLik -44189.9 edf 424.43 eps 0.1254 iteration   4
AICc 86452.65 logPost -2411216 logLik -42768.9 edf 415.34 eps 0.0628 iteration   5
AICc 86231.98 logPost -2411101 logLik -42661.9 edf 412.53 eps 0.0202 iteration   6
AICc 86215.73 logPost -2411097 logLik -42659.1 edf 408.13 eps 0.0027 iteration   7
AICc 86213.14 logPost -2411094 logLik -42658.2 edf 407.83 eps 0.0007 iteration   8
AICc 86211.64 logPost -2411092 logLik -42657.7 edf 407.62 eps 0.0003 iteration   9
AICc 86210.65 logPost -2411090 logLik -42657.4 edf 407.43 eps 0.0002 iteration  10
AICc 86205.44 logPost -2411208 logLik -42657.3 edf 405.37 eps 0.0001 iteration  11
AICc 86205.08 logPost -2411207 logLik -42657.3 edf 405.26 eps 0.0001 iteration  12
AICc 86202.45 logPost -2412410 logLik -42657.2 edf 404.20 eps 0.0001 iteration  13
AICc 86202.37 logPost -2412333 logLik -42657.2 edf 404.18 eps 0.0001 iteration  14
AICc 86202.37 logPost -2412333 logLik -42657.2 edf 404.18 eps 0.0001 iteration  14
elapsed time:  1.26min
Starting the sampler...

|                    |   0% 11.25min
|*                   |   5% 10.40min 32.84sec
|**                  |  10%  9.74min  1.08min
|***                 |  15%  8.80min  1.55min
|****                |  20%  8.47min  2.12min
|*****               |  25%  8.27min  2.76min
|******              |  30%  7.98min  3.42min
|*******             |  35%  7.58min  4.08min
|********            |  40%  7.09min  4.72min
|*********           |  45%  6.50min  5.32min
|**********          |  50%  6.00min  6.00min
|***********         |  55%  5.46min  6.67min
|************        |  60%  4.89min  7.33min
|*************       |  65%  4.30min  7.99min
|**************      |  70%  3.69min  8.62min
|***************     |  75%  3.09min  9.28min
|****************    |  80%  2.48min  9.93min
|*****************   |  85%  1.86min 10.56min
|******************  |  90%  1.24min 11.17min
|******************* |  95% 37.12sec 11.75min
|********************| 100%  0.00sec 12.37min
Code 
summary(multivariate_geo_sow_MRF_corr_endo)

Call:
bamlss(formula = f_rice_wheat_yield_MRF_corr_endo, family = mvnchol_bamlss(k = 4), 
    data = Irrig_Rev_rice_wheat, type = "modified")
---
Family: mvnchol 
Link function: mu1 = identity, mu2 = identity, mu3 = identity, mu4 = identity, lamdiag1 = log, lamdiag2 = log, lamdiag3 = log, lamdiag4 = log, lambda12 = identity, lambda13 = identity, lambda14 = identity, lambda23 = identity, lambda24 = identity, lambda34 = identity
*---
Formula mu1:
---
sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) + 
    s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + 
    s(sd_onset_82_15) + s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
                           Mean   2.5%    50%  97.5% parameters
(Intercept)              68.657 65.068 69.085 70.498      2.030
rice_duration_class_long -1.365 -2.006 -1.365 -0.714     -1.341
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(gw_2017).tau21            4.166e+00 1.290e-04 3.096e-02 3.818e+01    428.875
s(gw_2017).alpha            1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(gw_2017).edf              1.465e+00 9.987e-01 1.023e+00 3.953e+00      7.459
s(onset_2017).tau21         9.109e+00 1.164e-04 2.931e-01 8.568e+01   3370.676
s(onset_2017).alpha         1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(onset_2017).edf           1.809e+00 9.986e-01 1.298e+00 4.857e+00      8.245
s(monsoon_onset_dev).tau21  7.943e+00 2.625e-04 1.758e-01 7.388e+01     16.446
s(monsoon_onset_dev).alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(monsoon_onset_dev).edf    1.955e+00 9.997e-01 1.241e+00 5.247e+00      3.687
s(median_onset_82_15).tau21 2.851e+01 7.210e-04 8.607e+00 1.819e+02     18.779
s(median_onset_82_15).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(median_onset_82_15).edf   3.125e+00 1.001e+00 3.134e+00 6.203e+00      3.766
s(sd_onset_82_15).tau21     7.979e+00 8.517e-04 5.694e-01 6.914e+01      0.537
s(sd_onset_82_15).alpha     1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(sd_onset_82_15).edf       1.823e+00 1.001e+00 1.396e+00 4.743e+00      1.379
s(District,id='mrf1').tau21 3.093e+01 1.580e+01 2.913e+01 5.584e+01   1484.377
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='mrf1').edf   3.430e+01 3.394e+01 3.431e+01 3.457e+01     34.979
s(District,id='re2').tau21  1.606e+04 1.019e+04 1.538e+04 2.517e+04      1.087
s(District,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='re2').edf    3.599e+01 3.599e+01 3.599e+01 3.600e+01     35.240
---
Formula mu2:
---
b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + 
    s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + 
    s(p2o5perha_rice) + s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
                            Mean    2.5%     50%   97.5% parameters
(Intercept)              2.16338 0.43280 1.97697 3.90940      0.000
rice_duration_class_long 0.10365 0.03548 0.10337 0.17617      0.096
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                       Mean      2.5%       50%     97.5%
s(Res_rice_sow).tau21             2.988e+00 8.347e-02 1.431e+00 1.520e+01
s(Res_rice_sow).alpha             1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(Res_rice_sow).edf               5.567e+00 2.711e+00 5.607e+00 8.120e+00
s(g_q5305_irrig_times_rice).tau21 1.784e+00 2.339e-01 1.343e+00 6.251e+00
s(g_q5305_irrig_times_rice).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(g_q5305_irrig_times_rice).edf   6.203e+00 4.452e+00 6.275e+00 7.714e+00
s(nperha_rice).tau21              8.336e-01 6.327e-03 4.352e-01 4.123e+00
s(nperha_rice).alpha              1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(nperha_rice).edf                4.117e+00 1.326e+00 4.162e+00 6.756e+00
s(p2o5perha_rice).tau21           1.930e-01 1.197e-04 2.762e-02 1.220e+00
s(p2o5perha_rice).alpha           1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(p2o5perha_rice).edf             2.307e+00 1.009e+00 1.993e+00 5.123e+00
s(District,id='mrf1').tau21       3.021e-03 7.407e-05 9.868e-04 1.600e-02
s(District,id='mrf1').alpha       1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='mrf1').edf         1.459e+01 1.845e+00 1.295e+01 3.014e+01
s(District,id='re2').tau21        4.972e+00 1.205e-01 4.133e+00 1.458e+01
s(District,id='re2').alpha        1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='re2').edf          3.529e+01 3.319e+01 3.577e+01 3.593e+01
                                  parameters
s(Res_rice_sow).tau21                 16.244
s(Res_rice_sow).alpha                     NA
s(Res_rice_sow).edf                    8.170
s(g_q5305_irrig_times_rice).tau21      1.698
s(g_q5305_irrig_times_rice).alpha         NA
s(g_q5305_irrig_times_rice).edf        6.538
s(nperha_rice).tau21                   0.991
s(nperha_rice).alpha                      NA
s(nperha_rice).edf                     5.103
s(p2o5perha_rice).tau21                0.172
s(p2o5perha_rice).alpha                   NA
s(p2o5perha_rice).edf                  3.239
s(District,id='mrf1').tau21            0.035
s(District,id='mrf1').alpha               NA
s(District,id='mrf1').edf             32.259
s(District,id='re2').tau21             8.948
s(District,id='re2').alpha                NA
s(District,id='re2').edf              35.894
---
Formula mu3:
---
sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) + 
    s(District, bs = "mrf", xt = list(penalty = K)) + s(District, 
    bs = "re")
-
Parametric coefficients:
                     Mean    2.5%     50%   97.5% parameters
(Intercept)       114.653 109.482 114.477 120.585      3.830
variety_type_NMWV  -3.072  -3.739  -3.074  -2.457     -3.073
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(harvest_day_rice).tau21     228.447    44.418   160.009   703.064   1035.210
s(harvest_day_rice).alpha       1.000     1.000     1.000     1.000         NA
s(harvest_day_rice).edf         5.935     4.527     5.900     7.390      8.391
s(District,id='mrf1').tau21    53.987     9.438    45.562   154.096   4326.837
s(District,id='mrf1').alpha     1.000     1.000     1.000     1.000         NA
s(District,id='mrf1').edf      34.590    33.753    34.685    34.847     34.994
s(District,id='re2').tau21  49904.680 32149.456 48242.707 80335.112      1.087
s(District,id='re2').alpha      1.000     1.000     1.000     1.000         NA
s(District,id='re2').edf       35.999    35.998    35.999    35.999     35.240
---
Formula mu4:
---
l_ton_per_hectare ~ 1 + variety_type_NMWV + s(Res_wheat_sow) + 
    g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", 
    xt = list(penalty = K)) + s(District, bs = "re")
-
Parametric coefficients:
                          Mean       2.5%        50%      97.5% parameters
(Intercept)         -0.3451323 -0.9641137 -0.4141525  0.3765946     -1.419
variety_type_NMWV    0.3422674  0.2821928  0.3433422  0.3966434      0.337
g_q5305_irrig_times  0.4219312  0.3934013  0.4220813  0.4499296      0.424
nperha               0.0015495  0.0009506  0.0015567  0.0021201      0.002
p2o5perha            0.0025588  0.0014389  0.0025537  0.0036592      0.002
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(Res_wheat_sow).tau21      2.119e-02 7.832e-05 1.797e-03 1.665e-01      0.238
s(Res_wheat_sow).alpha      1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(Res_wheat_sow).edf        1.602e+00 1.011e+00 1.222e+00 4.012e+00      4.795
s(District,id='mrf1').tau21 1.002e-02 5.312e-05 5.624e-03 4.261e-02      0.010
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='mrf1').edf   2.528e+01 2.964e+00 2.937e+01 3.356e+01     31.282
s(District,id='re2').tau21  3.750e+00 1.161e+00 3.640e+00 7.587e+00      5.082
s(District,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='re2').edf    3.586e+01 3.568e+01 3.588e+01 3.594e+01     35.915
---
Formula lamdiag1:
---
lamdiag1 ~ 1
-
Parametric coefficients:
              Mean   2.5%    50%  97.5% parameters
(Intercept) -2.216 -2.236 -2.216 -2.194     -2.207
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9897 0.9057 1.0000     1
---
Formula lamdiag2:
---
lamdiag2 ~ 1
-
Parametric coefficients:
               Mean    2.5%     50%   97.5% parameters
(Intercept) 0.06659 0.04544 0.06608 0.08778      0.074
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9913 0.9201 1.0000     1
---
Formula lamdiag3:
---
lamdiag3 ~ 1
-
Parametric coefficients:
              Mean   2.5%    50%  97.5% parameters
(Intercept) -2.145 -2.167 -2.145 -2.126     -2.141
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9896 0.9165 1.0000     1
---
Formula lamdiag4:
---
lamdiag4 ~ 1
-
Parametric coefficients:
              Mean   2.5%    50%  97.5% parameters
(Intercept) 0.5177 0.4980 0.5172 0.5393      0.527
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9910 0.9193 1.0000     1
---
Formula lambda12:
---
lambda12 ~ 1
-
Parametric coefficients:
                Mean     2.5%      50%    97.5% parameters
(Intercept) -0.01094 -0.04337 -0.00753  0.02846     -0.004
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda13:
---
lambda13 ~ 1
-
Parametric coefficients:
                Mean     2.5%      50%    97.5% parameters
(Intercept) 0.007317 0.003840 0.007325 0.010682      0.008
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda14:
---
lambda14 ~ 1
-
Parametric coefficients:
                Mean     2.5%      50%    97.5% parameters
(Intercept) 0.011075 0.002135 0.010954 0.018862      0.011
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda23:
---
lambda23 ~ 1
-
Parametric coefficients:
                 Mean      2.5%       50%     97.5% parameters
(Intercept)  0.023666 -0.007654  0.023488  0.056214      0.027
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda24:
---
lambda24 ~ s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
               Mean    2.5%     50%   97.5% parameters
(Intercept) -0.2383 -0.2905 -0.2362 -0.1995     -0.036
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(District,id='mrf1').tau21 2.359e-03 9.909e-05 2.021e-03 6.741e-03      0.005
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='mrf1').edf   1.464e+01 1.841e+00 1.584e+01 2.434e+01     22.153
s(District,id='re2').tau21  9.752e-03 1.593e-04 6.349e-03 3.330e-02      0.024
s(District,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='re2').edf    1.392e+01 6.479e-01 1.434e+01 2.683e+01     24.454
---
Formula lambda34:
---
lambda34 ~ 1
-
Parametric coefficients:
                Mean     2.5%      50%    97.5% parameters
(Intercept) -0.07590 -0.11841 -0.08114 -0.01123     -0.009
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Sampler summary:
-
DIC = 85815.14 logLik = -42796.87 pd = 221.3947
runtime = 745.15
---
Optimizer summary:
-
AICc = 86202.37 edf = 404.1805 logLik = -42657.27
logPost = -2412333 nobs = 4527 runtime = 75.49
Code 
nd <- data.frame("District" = unique(Irrig_Rev_rice_wheat$District))

# Focusing on the cross-equation correlations

## Predict for the structured spatial effects.
p_str_multivariate_geo_sow_MRF_corr_endo_rice_y <- predict(multivariate_geo_sow_MRF_corr_endo, newdata = nd, term = "s(District,id='mrf1')", intercept = FALSE)

p_str_multivariate_geo_sow_MRF_corr_endo_rice_ydt <- as.data.frame(p_str_multivariate_geo_sow_MRF_corr_endo_rice_y)


## And the unstructured spatial effect.
p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y <- predict(multivariate_geo_sow_MRF_corr_endo, newdata = nd, term = "s(District,id='re2')", intercept = FALSE)

# Rice sowing spatial equation

plotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu1, id = nd$District, main = expression(mu(rice_sowing)), title = "Structured spatial effect")

Code 
plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu1, id = nd$District, main = expression(mu(rice_sowing)), title = "Unstructured spatial effect")

Code 
# Rice yield spatial equation
plotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu2, id = nd$District, main = expression(mu(rice_yield)), title = "Structured spatial effect")

Code 
plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu2, id = nd$District, main = expression(mu(rice_yield)), title = "Unstructured spatial effect")

Code 
# Wheat sowing equation
plotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu3, id = nd$District, main = expression(mu(wheat_sowing)), title = "Structured spatial effect")

Code 
plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu3, id = nd$District, main = expression(mu(wheat_sowing)), title = "Unstructured spatial effect")

Code 
# Wheat yield spatial equation
plotmap(India_aoi_sp_bnd, x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$mu4, id = nd$District, main = expression(mu(wheat_yield)), title = "Structured spatial effect")

Code 
plotmap(India_aoi_sp_bnd, x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$mu4, id = nd$District, main = expression(mu(wheat_yield)), title = "Unstructured spatial effect")

Code 
# Focusing on the cross-equation correlations

## MRF smooth effect.
plotmap(India_aoi_sp_bnd,
    x = p_str_multivariate_geo_sow_MRF_corr_endo_rice_y$lambda24, id = nd$District,
    main = expression(lambda(rice, wheat)), title = "Structured spatial effect"
)

Code 
## Random effects.
plotmap(India_aoi_sp_bnd,
    x = p_unstr_multivariate_geo_sow_MRF_corr_endo_rice_y$lambda24, id = nd$District, main = expression(lambda(rice, wheat)), title = "Unstructured spatial effect"
)

Code 
# Rice sowing equation : Non-linear relationships
# s(gw_2017) + s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + s(sd_onset_82_15)
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu1", term = "s(gw_2017)")

Code 
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu1", term = "s(onset_2017)")

Code 
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu1", term = "s(monsoon_onset_dev)")

Code 
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu1", term = "s(median_onset_82_15)")

Code 
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu1", term = "s(sd_onset_82_15)")

Code 
# Rice yield equation
# s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice)

plot(multivariate_geo_sow_MRF_corr_endo, model = "mu2", term = "s(Res_rice_sow)")

Code 
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu2", term = "s(g_q5305_irrig_times_rice)")

Code 
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu2", term = "s(nperha_rice)")

Code 
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu2", term = "s(p2o5perha_rice)")

Code 
# Wheat sowing equation
# s(harvest_day_rice) + s(gw_2018)
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu3", term = "s(harvest_day_rice)")

Code 
# (multivariate_geo_sow_MRF_corr_endo, model = "mu3", term = "s(gw_2018)")

# Wheat yield equation
plot(multivariate_geo_sow_MRF_corr_endo, model = "mu4", term = "s(Res_wheat_sow)")

Code 
# Fitted values
multivariate_geo_sow_MRF_corr_endo_fitted_values <- multivariate_geo_sow_MRF_corr_endo$fitted

multivariate_geo_sow_MRF_corr_endo_fitted_values <- as.data.frame(multivariate_geo_sow_MRF_corr_endo_fitted_values)

# Merge the fitted results to the data and export



Irrig_Rev_rice_wheat_Mult_Res <- cbind(Irrig_Rev_rice_wheat, multivariate_geo_sow_MRF_corr_endo_fitted_values)


summary(lm(sowdate_fmt_rice_day ~ mu1, Irrig_Rev_rice_wheat_Mult_Res))

Call:
lm(formula = sowdate_fmt_rice_day ~ mu1, data = Irrig_Rev_rice_wheat_Mult_Res)

Residuals:
    Min      1Q  Median      3Q     Max 
-44.574  -5.397  -0.269   4.986  44.985 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.03383    3.41385   -0.01    0.992    
mu1          1.00018    0.01771   56.48   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.093 on 4525 degrees of freedom
Multiple R-squared:  0.4135,    Adjusted R-squared:  0.4134 
F-statistic:  3190 on 1 and 4525 DF,  p-value: < 2.2e-16
Code 
summary(lm(b_grain_yield_ton_per_ha_rice ~ l_ton_per_hectare, Irrig_Rev_rice_wheat_Mult_Res))

Call:
lm(formula = b_grain_yield_ton_per_ha_rice ~ l_ton_per_hectare, 
    data = Irrig_Rev_rice_wheat_Mult_Res)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2891 -0.7976 -0.0052  0.7454  3.3894 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        3.34109    0.05974   55.93   <2e-16 ***
l_ton_per_hectare  0.22766    0.01958   11.63   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.08 on 4525 degrees of freedom
Multiple R-squared:  0.02901,   Adjusted R-squared:  0.0288 
F-statistic: 135.2 on 1 and 4525 DF,  p-value: < 2.2e-16
Code 
plot((lm(sowdate_fmt_rice_day ~ mu1, Irrig_Rev_rice_wheat_Mult_Res)))

Code 
Irrig_Rev_rice_wheat_Mult_Res_sp <- Irrig_Rev_rice_wheat_Mult_Res
coordinates(Irrig_Rev_rice_wheat_Mult_Res_sp) <- c("o_largest_plot_gps_longitude", "o_largest_plot_gps_latitude")

# Map the correlations
# library(modelsummary)
# Irrig_Rev_rice_wheat_Mult_Res_dist <- datasummary(Heading("District") * District ~ Heading("N obs") * N + Heading("%") * Percent() + lambda24 * (Mean + SD), data = Irrig_Rev_rice_wheat_Mult_Res, output = "data.frame")

# library(dplyr)
# Irrig_Rev_rice_wheat_Mult_Res_dist <- rename(Irrig_Rev_rice_wheat_Mult_Res_dist, "Mean_Rice_Wheat_Rho" = "Mean")
# Irrig_Rev_rice_wheat_Mult_Res_dist <- rename(Irrig_Rev_rice_wheat_Mult_Res_dist, "SD_Rice_Wheat_Rho" = "SD")

# Irrig_Rev_rice_wheat_Mult_Res_dist <- subset(Irrig_Rev_rice_wheat_Mult_Res_dist, Irrig_Rev_rice_wheat_Mult_Res_dist$District != "Purnia")

# rice_wheat_yield_rho_dist_sf <- merge(India_aoi_sf, Irrig_Rev_rice_wheat_Mult_Res_dist, by = "District")

# rice_wheat_yield_rho_dist_sf$Mean_Rice_Wheat_Rho <- as.numeric(rice_wheat_yield_rho_dist_sf$Mean_Rice_Wheat_Rho)
# library(mapview)
# mapview(rice_wheat_yield_rho_dist_sf, zcol = "Mean_Rice_Wheat_Rho", layer.name = "Rice wheat equation correlation")
# library(sf)
# rice_wheat_yield_rho_dist_sf_sp <- as_Spatial(rice_wheat_yield_rho_dist_sf)
# library(tmap)
# tmap_mode("view")
# rice_wheat_yield_rho_dist_sf_sp_map <- tm_shape(rice_wheat_yield_rho_dist_sf_sp) +
#     tm_polygons(col = "Mean_Rice_Wheat_Rho", title = "Rice wheat equation correlation", style = "quantile") +
#     tm_layout(legend.outside = TRUE)

# tmap_save(rice_wheat_yield_rho_dist_sf_sp_map, "figures/rice_wheat_yield_rho_dist_sf_sp_map .png")

5 Factors affecting the correlation structure

Code 
f_rice_wheat_yield_MRF_corr_endo_General <- list(
    sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) + s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + s(sd_onset_82_15) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    l_ton_per_hectare ~ 1 + variety_type_NMWV + s(Res_wheat_sow) + g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", xt = list("penalty" = K)) + s(District, bs = "re"),
    lamdiag1 ~ 1,
    lamdiag2 ~ 1,
    lamdiag3 ~ 1,
    lamdiag4 ~ 1,
    lambda12 ~ 1,
    lambda13 ~ 1,
    lambda14 ~ 1,
    lambda23 ~ 1,
    lambda24 ~ s(Res_rice_sow)+s(Res_wheat_sow)+rice_duration_class_long + s(gw_2017) + s(onset_2017)+variety_type_NMWV+s(District, bs = "mrf", xt = list("penalty" = K)) + s(District, bs = "re"),
    lambda34 ~ 1
)

multivariate_geo_sow_MRF_corr_endo_General <- bamlss(f_rice_wheat_yield_MRF_corr_endo_General, type = "modified", family = mvnchol_bamlss(k = 4), data = Irrig_Rev_rice_wheat)
AICc 349496.0 logPost -2543650 logLik -174245. edf 451.97 eps 1.0000 iteration   1
AICc 157756.9 logPost -2447315 logLik -78369.9 edf 457.10 eps 0.2663 iteration   2
AICc 104296.5 logPost -2420163 logLik -51655.0 edf 444.64 eps 0.2173 iteration   3
AICc 89324.62 logPost -2412631 logLik -44177.5 edf 437.79 eps 0.1257 iteration   4
AICc 86460.20 logPost -2411212 logLik -42756.4 edf 428.65 eps 0.0632 iteration   5
AICc 86239.60 logPost -2411096 logLik -42649.5 edf 425.88 eps 0.0203 iteration   6
AICc 86223.26 logPost -2411092 logLik -42646.7 edf 421.47 eps 0.0027 iteration   7
AICc 86220.64 logPost -2411088 logLik -42645.8 edf 421.15 eps 0.0008 iteration   8
AICc 86219.11 logPost -2411086 logLik -42645.3 edf 420.93 eps 0.0004 iteration   9
AICc 86218.10 logPost -2411084 logLik -42645.0 edf 420.73 eps 0.0002 iteration  10
AICc 86212.87 logPost -2411201 logLik -42644.9 edf 418.67 eps 0.0002 iteration  11
AICc 86210.00 logPost -2412487 logLik -42644.9 edf 417.53 eps 0.0002 iteration  12
AICc 86209.82 logPost -2412474 logLik -42644.9 edf 417.48 eps 0.0001 iteration  13
AICc 86209.73 logPost -2412389 logLik -42644.8 edf 417.45 eps 0.0001 iteration  14
AICc 86209.64 logPost -2412389 logLik -42644.8 edf 417.42 eps 0.0001 iteration  15
AICc 86209.64 logPost -2412389 logLik -42644.8 edf 417.42 eps 0.0001 iteration  15
elapsed time:  1.67min
Starting the sampler...

|                    |   0% 11.50min
|*                   |   5% 11.36min 35.86sec
|**                  |  10% 10.54min  1.17min
|***                 |  15%  9.83min  1.74min
|****                |  20%  9.63min  2.41min
|*****               |  25%  9.46min  3.15min
|******              |  30%  9.02min  3.87min
|*******             |  35%  8.52min  4.59min
|********            |  40%  7.93min  5.29min
|*********           |  45%  7.32min  5.99min
|**********          |  50%  6.69min  6.69min
|***********         |  55%  6.05min  7.39min
|************        |  60%  5.40min  8.09min
|*************       |  65%  4.74min  8.80min
|**************      |  70%  4.07min  9.50min
|***************     |  75%  3.40min 10.20min
|****************    |  80%  2.73min 10.91min
|*****************   |  85%  2.06min 11.67min
|******************  |  90%  1.38min 12.45min
|******************* |  95% 41.67sec 13.20min
|********************| 100%  0.00sec 13.91min
Code 
summary(multivariate_geo_sow_MRF_corr_endo_General)

Call:
bamlss(formula = f_rice_wheat_yield_MRF_corr_endo_General, family = mvnchol_bamlss(k = 4), 
    data = Irrig_Rev_rice_wheat, type = "modified")
---
Family: mvnchol 
Link function: mu1 = identity, mu2 = identity, mu3 = identity, mu4 = identity, lamdiag1 = log, lamdiag2 = log, lamdiag3 = log, lamdiag4 = log, lambda12 = identity, lambda13 = identity, lambda14 = identity, lambda23 = identity, lambda24 = identity, lambda34 = identity
*---
Formula mu1:
---
sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) + 
    s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + 
    s(sd_onset_82_15) + s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
                            Mean    2.5%     50%   97.5% parameters
(Intercept)              69.2044 66.0664 69.0422 73.0918      2.030
rice_duration_class_long -1.3425 -1.9950 -1.3446 -0.6804     -1.339
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(gw_2017).tau21            1.384e+01 9.369e-05 4.354e-01 1.444e+02    422.419
s(gw_2017).alpha            1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(gw_2017).edf              1.893e+00 9.981e-01 1.258e+00 5.384e+00      7.444
s(onset_2017).tau21         3.598e+00 2.236e-04 1.113e-01 2.751e+01   3149.404
s(onset_2017).alpha         1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(onset_2017).edf           1.497e+00 9.995e-01 1.127e+00 3.753e+00      8.208
s(monsoon_onset_dev).tau21  1.420e+01 1.552e-03 2.525e+00 9.327e+01     15.892
s(monsoon_onset_dev).alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(monsoon_onset_dev).edf    2.534e+00 1.002e+00 2.259e+00 5.423e+00      3.656
s(median_onset_82_15).tau21 4.901e+00 2.852e-04 1.789e-02 5.389e+01     18.463
s(median_onset_82_15).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(median_onset_82_15).edf   1.546e+00 9.998e-01 1.024e+00 4.762e+00      3.751
s(sd_onset_82_15).tau21     2.742e+00 1.451e-04 3.032e-02 2.532e+01      0.501
s(sd_onset_82_15).alpha     1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(sd_onset_82_15).edf       1.410e+00 9.989e-01 1.029e+00 3.681e+00      1.359
s(District,id='mrf1').tau21 1.848e+01 6.094e+00 1.647e+01 4.042e+01   1473.762
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='mrf1').edf   3.389e+01 3.303e+01 3.396e+01 3.445e+01     34.979
s(District,id='re2').tau21  1.636e+04 1.023e+04 1.566e+04 2.643e+04      1.087
s(District,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='re2').edf    3.599e+01 3.599e+01 3.599e+01 3.600e+01     35.240
---
Formula mu2:
---
b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + 
    s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + 
    s(p2o5perha_rice) + s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
                            Mean    2.5%     50%   97.5% parameters
(Intercept)              1.25217 0.43813 1.14065 2.42250     -0.001
rice_duration_class_long 0.10039 0.03315 0.09851 0.16726      0.098
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                       Mean      2.5%       50%     97.5%
s(Res_rice_sow).tau21             3.436e+00 9.838e-02 1.475e+00 1.706e+01
s(Res_rice_sow).alpha             1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(Res_rice_sow).edf               5.675e+00 2.846e+00 5.701e+00 8.189e+00
s(g_q5305_irrig_times_rice).tau21 1.800e+00 2.744e-01 1.290e+00 6.920e+00
s(g_q5305_irrig_times_rice).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(g_q5305_irrig_times_rice).edf   6.226e+00 4.614e+00 6.233e+00 7.781e+00
s(nperha_rice).tau21              6.891e-01 1.623e-02 3.422e-01 3.336e+00
s(nperha_rice).alpha              1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(nperha_rice).edf                3.946e+00 1.674e+00 3.921e+00 6.527e+00
s(p2o5perha_rice).tau21           1.232e-01 1.389e-04 6.377e-03 1.218e+00
s(p2o5perha_rice).alpha           1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(p2o5perha_rice).edf             1.927e+00 1.010e+00 1.370e+00 5.094e+00
s(District,id='mrf1').tau21       1.366e-03 4.887e-05 6.491e-04 6.166e-03
s(District,id='mrf1').alpha       1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='mrf1').edf         1.155e+01 1.279e+00 1.008e+01 2.573e+01
s(District,id='re2').tau21        8.281e+00 2.283e+00 7.833e+00 1.717e+01
s(District,id='re2').alpha        1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='re2').edf          3.585e+01 3.562e+01 3.587e+01 3.594e+01
                                  parameters
s(Res_rice_sow).tau21                 17.005
s(Res_rice_sow).alpha                     NA
s(Res_rice_sow).edf                    8.216
s(g_q5305_irrig_times_rice).tau21      1.765
s(g_q5305_irrig_times_rice).alpha         NA
s(g_q5305_irrig_times_rice).edf        6.581
s(nperha_rice).tau21                   1.007
s(nperha_rice).alpha                      NA
s(nperha_rice).edf                     5.126
s(p2o5perha_rice).tau21                0.180
s(p2o5perha_rice).alpha                   NA
s(p2o5perha_rice).edf                  3.292
s(District,id='mrf1').tau21            0.035
s(District,id='mrf1').alpha               NA
s(District,id='mrf1').edf             32.235
s(District,id='re2').tau21             8.974
s(District,id='re2').alpha                NA
s(District,id='re2').edf              35.893
---
Formula mu3:
---
sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) + 
    s(District, bs = "mrf", xt = list(penalty = K)) + s(District, 
    bs = "re")
-
Parametric coefficients:
                     Mean    2.5%     50%   97.5% parameters
(Intercept)       110.043 107.192 110.123 112.566      3.830
variety_type_NMWV  -3.067  -3.716  -3.071  -2.442     -3.073
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(harvest_day_rice).tau21     232.822    45.619   165.884   752.545   1027.121
s(harvest_day_rice).alpha       1.000     1.000     1.000     1.000         NA
s(harvest_day_rice).edf         5.950     4.434     5.907     7.447      8.387
s(District,id='mrf1').tau21   168.059    80.428   152.662   327.771   4298.849
s(District,id='mrf1').alpha     1.000     1.000     1.000     1.000         NA
s(District,id='mrf1').edf      34.867    34.741    34.870    34.961     34.994
s(District,id='re2').tau21  52604.006 33225.833 50703.021 84863.013      1.087
s(District,id='re2').alpha      1.000     1.000     1.000     1.000         NA
s(District,id='re2').edf       35.999    35.998    35.999    35.999     35.240
---
Formula mu4:
---
l_ton_per_hectare ~ 1 + variety_type_NMWV + s(Res_wheat_sow) + 
    g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", 
    xt = list(penalty = K)) + s(District, bs = "re")
-
Parametric coefficients:
                          Mean       2.5%        50%      97.5% parameters
(Intercept)         -0.4432149 -1.1226803 -0.4751837  0.5361674     -1.417
variety_type_NMWV    0.3438003  0.2935946  0.3435507  0.3932263      0.337
g_q5305_irrig_times  0.4227435  0.3952499  0.4232928  0.4494365      0.425
nperha               0.0015125  0.0009256  0.0015230  0.0020544      0.001
p2o5perha            0.0025719  0.0014937  0.0025608  0.0037184      0.003
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(Res_wheat_sow).tau21      1.091e-02 7.030e-05 1.574e-03 1.004e-01      0.179
s(Res_wheat_sow).alpha      1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(Res_wheat_sow).edf        1.453e+00 1.010e+00 1.192e+00 3.490e+00      4.475
s(District,id='mrf1').tau21 1.130e-02 3.038e-04 7.371e-03 3.612e-02      0.010
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='mrf1').edf   2.799e+01 1.060e+01 3.038e+01 3.341e+01     31.273
s(District,id='re2').tau21  4.202e+00 9.197e-01 3.924e+00 9.293e+00      5.097
s(District,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='re2').edf    3.586e+01 3.561e+01 3.589e+01 3.595e+01     35.915
---
Formula lamdiag1:
---
lamdiag1 ~ 1
-
Parametric coefficients:
              Mean   2.5%    50%  97.5% parameters
(Intercept) -2.216 -2.236 -2.216 -2.196     -2.207
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9893 0.8980 1.0000     1
---
Formula lamdiag2:
---
lamdiag2 ~ 1
-
Parametric coefficients:
               Mean    2.5%     50%   97.5% parameters
(Intercept) 0.06608 0.04506 0.06609 0.08676      0.074
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9912 0.9291 1.0000     1
---
Formula lamdiag3:
---
lamdiag3 ~ 1
-
Parametric coefficients:
              Mean   2.5%    50%  97.5% parameters
(Intercept) -2.146 -2.165 -2.146 -2.123     -2.141
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9892 0.9017 0.9999     1
---
Formula lamdiag4:
---
lamdiag4 ~ 1
-
Parametric coefficients:
              Mean   2.5%    50%  97.5% parameters
(Intercept) 0.5192 0.4982 0.5194 0.5397      0.529
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9910 0.9196 1.0000     1
---
Formula lambda12:
---
lambda12 ~ 1
-
Parametric coefficients:
                  Mean       2.5%        50%      97.5% parameters
(Intercept)  0.0011779 -0.0188817  0.0009202  0.0275793     -0.004
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda13:
---
lambda13 ~ 1
-
Parametric coefficients:
                Mean     2.5%      50%    97.5% parameters
(Intercept) 0.007509 0.003970 0.007520 0.011361      0.008
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda14:
---
lambda14 ~ 1
-
Parametric coefficients:
                Mean     2.5%      50%    97.5% parameters
(Intercept) 0.008838 0.003256 0.008796 0.013986      0.011
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda23:
---
lambda23 ~ 1
-
Parametric coefficients:
                 Mean      2.5%       50%     97.5% parameters
(Intercept)  0.023089 -0.009292  0.023265  0.056951      0.027
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda24:
---
lambda24 ~ s(Res_rice_sow) + s(Res_wheat_sow) + rice_duration_class_long + 
    s(gw_2017) + s(onset_2017) + variety_type_NMWV + s(District, 
    bs = "mrf", xt = list(penalty = K)) + s(District, bs = "re")
-
Parametric coefficients:
                             Mean     2.5%      50%    97.5% parameters
(Intercept)              -0.17848 -0.24809 -0.17781 -0.10775      0.005
rice_duration_class_long -0.02497 -0.10360 -0.02418  0.04930     -0.019
variety_type_NMWV        -0.09328 -0.17084 -0.09252 -0.01647     -0.067
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(Res_rice_sow).tau21       1.529e-01 1.272e-04 1.007e-02 1.519e+00      1.693
s(Res_rice_sow).alpha       1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(Res_rice_sow).edf         1.803e+00 1.007e+00 1.395e+00 5.469e+00      5.426
s(Res_wheat_sow).tau21      1.700e-02 9.987e-05 2.104e-03 1.790e-01      0.012
s(Res_wheat_sow).alpha      1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(Res_wheat_sow).edf        1.304e+00 1.004e+00 1.088e+00 2.991e+00      1.393
s(gw_2017).tau21            5.304e-02 9.734e-05 5.875e-03 4.649e-01      0.001
s(gw_2017).alpha            1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(gw_2017).edf              1.599e+00 1.005e+00 1.265e+00 3.782e+00      1.042
s(onset_2017).tau21         2.070e-02 6.115e-05 1.994e-03 1.472e-01      1.275
s(onset_2017).alpha         1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(onset_2017).edf           1.400e+00 1.006e+00 1.171e+00 2.983e+00      4.996
s(District,id='mrf1').tau21 2.852e-03 8.583e-05 2.535e-03 7.671e-03      0.004
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='mrf1').edf   1.553e+01 1.557e+00 1.723e+01 2.486e+01     21.116
s(District,id='re2').tau21  7.156e-03 1.261e-04 3.300e-03 3.121e-02      0.023
s(District,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='re2').edf    1.092e+01 4.917e-01 9.199e+00 2.639e+01     24.188
---
Formula lambda34:
---
lambda34 ~ 1
-
Parametric coefficients:
                Mean     2.5%      50%    97.5% parameters
(Intercept) -0.06752 -0.13892 -0.07741  0.03476     -0.009
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Sampler summary:
-
DIC = 85815.34 logLik = -42794.18 pd = 226.9757
runtime = 837.55
---
Optimizer summary:
-
AICc = 86209.64 edf = 417.425 logLik = -42644.88
logPost = -2412389 nobs = 4527 runtime = 100.41
Code 
plot(multivariate_geo_sow_MRF_corr_endo_General, model = "lambda24", term = "s(gw_2017)")

Code 
plot(multivariate_geo_sow_MRF_corr_endo_General, model = "lambda24", term = "s(onset_2017)")

Code 
plot(multivariate_geo_sow_MRF_corr_endo_General, model = "lambda24", term = "s(Res_rice_sow)")

Code 
plot(multivariate_geo_sow_MRF_corr_endo_General, model = "lambda24", term = "s(Res_wheat_sow)")

Code 
#save.image("RW_System_SDS.RData")

6 Including Soils and other variables [For Bihar Only]

Code 
library(sp)
library(mgcv)
library(bamlss)

# Multivariate geoadditive model
# remotes::install_git("https://git.uibk.ac.at/c4031039/mvnchol")
# library(mvnchol)
library(BayesX)
library(R2BayesX)
library(sf)
library(spdep)
library(rio)

Irrig_Rev_rice_wheat <- import("Irrig_Rev_rice_wheat_Updated.csv")


Irrig_Rev_rice_wheat$harvest_day_rice <- Irrig_Rev_rice_wheat$l_crop_duration_days_rice + Irrig_Rev_rice_wheat$sowdate_fmt_rice_day

# Irrig_Rev_rice_wheat$harvest_day_wheat <- Irrig_Rev_rice_wheat$l_crop_duration_days + Irrig_Rev_rice_wheat$sowdate_fmt_wheat_day

# Remove NAs in the monsoon variables
Irrig_Rev_rice_wheat <- subset(Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$onset_2017)))
Irrig_Rev_rice_wheat <- subset(Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$monsoon_onset_dev)))
Irrig_Rev_rice_wheat <- subset(Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$median_onset_82_15)))
Irrig_Rev_rice_wheat <- subset(Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$sd_onset_82_15)))

Irrig_Rev_rice_wheat <- subset(Irrig_Rev_rice_wheat, !(is.na(Irrig_Rev_rice_wheat$May_tmax_17)))

shpname <- file.path(getwd(), "shp", "India_aoi_sf_sp")
India_aoi_sp_bnd <- BayesX::shp2bnd(shpname = shpname, regionnames = "District", check.is.in = F)
Reading map ... finished
Note: map consists originally of 50 polygons
Note: map consists of 47 regions
Code 
K <- neighbormatrix(India_aoi_sp_bnd)
head(K)
           Araria Arwal Aurangabad Banka Begusarai Bhagalpur Arah Buxar
Araria          4     0          0     0         0         0    0     0
Arwal           0     6         -1     0         0         0   -1     0
Aurangabad      0    -1          3     0         0         0    0     0
Banka           0     0          0     3         0        -1    0     0
Begusarai       0     0          0     0         5         0    0     0
Bhagalpur       0     0          0    -1         0         6    0     0
           Darbhanga Gaya Gopalganj Jamui Jehanabad Kaimur Katihar Khagaria
Araria             0    0         0     0         0      0       0        0
Arwal              0   -1         0     0        -1      0       0        0
Aurangabad         0   -1         0     0         0      0       0        0
Banka              0    0         0    -1         0      0       0        0
Begusarai          0    0         0     0         0      0       0       -1
Bhagalpur          0    0         0     0         0      0      -1       -1
           Kishanganj Lakhisarai Madhepura Madhubani Munger Muzaffarpur Nalanda
Araria             -1          0        -1         0      0           0       0
Arwal               0          0         0         0      0           0       0
Aurangabad          0          0         0         0      0           0       0
Banka               0          0         0         0     -1           0       0
Begusarai           0         -1         0         0     -1           0       0
Bhagalpur           0          0        -1         0     -1           0       0
           Nawada WestChamparan Patna EastChamparan Purnia Rohtas Saharsa
Araria          0             0     0             0     -1      0       0
Arwal           0             0    -1             0      0     -1       0
Aurangabad      0             0     0             0      0     -1       0
Banka           0             0     0             0      0      0       0
Begusarai       0             0    -1             0      0      0       0
Bhagalpur       0             0     0             0     -1      0       0
           Samastipur Saran Sheikhpura Sheohar Sitamarhi Siwan Supaul Vaishali
Araria              0     0          0       0         0     0     -1        0
Arwal               0     0          0       0         0     0      0        0
Aurangabad          0     0          0       0         0     0      0        0
Banka               0     0          0       0         0     0      0        0
Begusarai          -1     0          0       0         0     0      0        0
Bhagalpur           0     0          0       0         0     0      0        0
           Balia Chandauli Deoria Gazipur Gorakhpur Kushinagar Maharajganj Mau
Araria         0         0      0       0         0          0           0   0
Arwal          0         0      0       0         0          0           0   0
Aurangabad     0         0      0       0         0          0           0   0
Banka          0         0      0       0         0          0           0   0
Begusarai      0         0      0       0         0          0           0   0
Bhagalpur      0         0      0       0         0          0           0   0
           Siddharthnagar
Araria                  0
Arwal                   0
Aurangabad              0
Banka                   0
Begusarai               0
Bhagalpur               0
Code 
## Also need to transform to factor for
## setting up the MRF smooth.
Irrig_Rev_rice_wheat$District <- as.factor(Irrig_Rev_rice_wheat$a_q103_district)

## Now note that not all regions are observed,
## therefore we need to remove those regions
## from the penalty matrix
rn <- rownames(K)
lv <- levels(Irrig_Rev_rice_wheat$District)
i <- rn %in% lv
K <- K[i, i]

# First stage analytics
library(sp)
library(mgcv)
library(bamlss)

# Multivariate geoadditive model
# remotes::install_git("https://git.uibk.ac.at/c4031039/mvnchol")
# install_github("https://github.com/meteosimon/mvnchol")
library(mvnchol)


# Rice first stage
f_sow_rice_1st_stage <- list(
    sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) + s(May_tmax_17) + s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + s(sd_onset_82_15) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re")
)

f_sow_rice_1st_stage_MRF <- bamlss(f_sow_rice_1st_stage, data = Irrig_Rev_rice_wheat, family = "gaussian")
AICc 38052.33 logPost -164341. logLik -18964.8 edf 60.445 eps 0.1196 iteration   1
AICc 34024.10 logPost -30660.3 logLik -16928.3 edf 82.182 eps 0.0978 iteration   2
AICc 32841.19 logPost -18077.9 logLik -16320.4 edf 97.892 eps 0.0612 iteration   3
AICc 32589.71 logPost -16803.5 logLik -16196.6 edf 96.050 eps 0.0271 iteration   4
AICc 32524.67 logPost -16754.2 logLik -16166.4 edf 93.812 eps 0.0057 iteration   5
AICc 32489.66 logPost -16715.1 logLik -16151.5 edf 91.392 eps 0.0015 iteration   6
AICc 32472.35 logPost -16693.1 logLik -16144.4 edf 89.848 eps 0.0007 iteration   7
AICc 32462.95 logPost -16695.5 logLik -16140.5 edf 89.058 eps 0.0004 iteration   8
AICc 32457.62 logPost -16618.9 logLik -16138.2 edf 88.747 eps 0.0003 iteration   9
AICc 32454.30 logPost -16613.1 logLik -16136.8 edf 88.459 eps 0.0002 iteration  10
AICc 32451.84 logPost -16793.1 logLik -16136.0 edf 88.067 eps 0.0001 iteration  11
AICc 32450.40 logPost -16832.9 logLik -16135.5 edf 87.887 eps 0.0001 iteration  12
AICc 32450.40 logPost -16832.9 logLik -16135.5 edf 87.887 eps 0.0001 iteration  12
elapsed time:  5.67sec
Starting the sampler...

|                    |   0%  1.05min
|*                   |   5%  1.09min  3.44sec
|**                  |  10%  1.03min  6.86sec
|***                 |  15% 58.20sec 10.27sec
|****                |  20% 55.24sec 13.81sec
|*****               |  25% 53.07sec 17.69sec
|******              |  30% 50.73sec 21.74sec
|*******             |  35% 47.62sec 25.64sec
|********            |  40% 44.32sec 29.55sec
|*********           |  45% 40.66sec 33.27sec
|**********          |  50% 36.50sec 36.50sec
|***********         |  55% 32.96sec 40.28sec
|************        |  60% 29.68sec 44.52sec
|*************       |  65% 26.12sec 48.50sec
|**************      |  70% 22.64sec 52.83sec
|***************     |  75% 19.03sec 57.09sec
|****************    |  80% 15.32sec  1.02min
|*****************   |  85% 11.57sec  1.09min
|******************  |  90%  7.72sec  1.16min
|******************* |  95%  3.87sec  1.23min
|********************| 100%  0.00sec  1.29min
Code 
fitted_f_sow_rice_1st_stage_MRF <- f_sow_rice_1st_stage_MRF$fitted

Irrig_Rev_rice_wheat$Res_rice_sow <- Irrig_Rev_rice_wheat$sowdate_fmt_rice_day - fitted_f_sow_rice_1st_stage_MRF$mu

# Wheat first stage
f_sow_wheat_1st_stage <- list(
    sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) + s(gw_2017) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re")
)

f_sow_wheat_1st_stage_MRF <- bamlss(f_sow_wheat_1st_stage, data = Irrig_Rev_rice_wheat, family = "gaussian")
AICc 47969.70 logPost -480141. logLik -23935.1 edf 49.173 eps 0.2531 iteration   1
AICc 43699.60 logPost -63125.9 logLik -21823.2 edf 26.412 eps 0.0752 iteration   2
AICc 39383.89 logPost -23961.6 logLik -19623.6 edf 67.239 eps 0.0828 iteration   3
AICc 35746.92 logPost -18513.4 logLik -17796.6 edf 75.537 eps 0.0792 iteration   4
AICc 33169.34 logPost -16984.5 logLik -16505.5 edf 77.723 eps 0.0792 iteration   5
AICc 32057.88 logPost -16496.1 logLik -15948.9 edf 78.558 eps 0.0635 iteration   6
AICc 31891.59 logPost -17061.8 logLik -15865.5 edf 78.787 eps 0.0289 iteration   7
AICc 31887.96 logPost -17061.5 logLik -15863.4 edf 79.094 eps 0.0043 iteration   8
AICc 31887.78 logPost -17061.7 logLik -15863.2 edf 79.153 eps 0.0001 iteration   9
AICc 31887.78 logPost -17061.7 logLik -15863.2 edf 79.153 eps 0.0001 iteration   9
elapsed time:  2.19sec
Starting the sampler...

|                    |   0% 44.03sec
|*                   |   5% 44.46sec  2.34sec
|**                  |  10% 42.12sec  4.68sec
|***                 |  15% 40.01sec  7.06sec
|****                |  20% 38.12sec  9.53sec
|*****               |  25% 37.20sec 12.40sec
|******              |  30% 35.58sec 15.25sec
|*******             |  35% 33.54sec 18.06sec
|********            |  40% 31.35sec 20.90sec
|*********           |  45% 28.64sec 23.43sec
|**********          |  50% 26.25sec 26.25sec
|***********         |  55% 23.52sec 28.75sec
|************        |  60% 20.90sec 31.35sec
|*************       |  65% 18.27sec 33.93sec
|**************      |  70% 15.75sec 36.75sec
|***************     |  75% 13.11sec 39.32sec
|****************    |  80% 10.56sec 42.23sec
|*****************   |  85%  7.91sec 44.84sec
|******************  |  90%  5.27sec 47.45sec
|******************* |  95%  2.63sec 50.01sec
|********************| 100%  0.00sec 53.01sec
Code 
fitted_f_sow_wheat_1st_stage_MRF <- f_sow_wheat_1st_stage_MRF$fitted

Irrig_Rev_rice_wheat$Res_wheat_sow <- Irrig_Rev_rice_wheat$sowdate_fmt_wheat_day - fitted_f_sow_wheat_1st_stage_MRF$mu


# Second stage analytics
f_rice_wheat_yield_MRF_corr_endo_General <- list(
    sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) + s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + s(sd_onset_82_15) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + s(p2o5perha_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) + s(District, bs = "mrf", xt = list("penalty" = K)) +
        s(District, bs = "re"),
    l_ton_per_hectare ~ 1 + variety_type_NMWV + s(Res_wheat_sow) + g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", xt = list("penalty" = K)) + s(District, bs = "re"),
    lamdiag1 ~ 1,
    lamdiag2 ~ 1,
    lamdiag3 ~ 1,
    lamdiag4 ~ 1,
    lambda12 ~ 1,
    lambda13 ~ 1,
    lambda14 ~ 1,
    lambda23 ~ 1,
    lambda24 ~ s(Res_rice_sow) + s(Res_wheat_sow) + rice_duration_class_long + s(gw_2017) + s(onset_2017) + variety_type_NMWV + s(District, bs = "mrf", xt = list("penalty" = K)) + s(District, bs = "re"),
    lambda34 ~ 1
)

multivariate_geo_sow_MRF_corr_endo_General <- bamlss(f_rice_wheat_yield_MRF_corr_endo_General, type = "modified", family = mvnchol_bamlss(k = 4), data = Irrig_Rev_rice_wheat)
AICc 356623.5 logPost -2493100 logLik -177815. edf 446.07 eps 1.0000 iteration   1
AICc 155195.3 logPost -2392042 logLik -77098.1 edf 449.01 eps 0.2813 iteration   2
AICc 102570.2 logPost -2365223 logLik -50795.9 edf 440.59 eps 0.2027 iteration   3
AICc 87819.30 logPost -2357820 logLik -43429.0 edf 433.66 eps 0.1178 iteration   4
AICc 84994.68 logPost -2356458 logLik -42027.1 edf 425.13 eps 0.0622 iteration   5
AICc 84777.36 logPost -2356772 logLik -41921.8 edf 422.38 eps 0.0204 iteration   6
AICc 84758.08 logPost -2356768 logLik -41919.3 edf 416.52 eps 0.0023 iteration   7
AICc 84756.01 logPost -2356765 logLik -41918.6 edf 416.29 eps 0.0007 iteration   8
AICc 84754.22 logPost -2356763 logLik -41918.2 edf 415.88 eps 0.0004 iteration   9
AICc 84753.43 logPost -2356762 logLik -41918.0 edf 415.72 eps 0.0002 iteration  10
AICc 84748.01 logPost -2356884 logLik -41917.9 edf 413.58 eps 0.0002 iteration  11
AICc 84747.69 logPost -2361228 logLik -41917.8 edf 413.48 eps 0.0002 iteration  12
AICc 84747.52 logPost -2361227 logLik -41917.8 edf 413.43 eps 0.0002 iteration  13
AICc 84747.39 logPost -2361227 logLik -41917.8 edf 413.39 eps 0.0002 iteration  14
AICc 84747.30 logPost -2361227 logLik -41917.8 edf 413.36 eps 0.0003 iteration  15
AICc 84747.24 logPost -2361226 logLik -41917.8 edf 413.33 eps 0.0002 iteration  16
AICc 84747.19 logPost -2361226 logLik -41917.8 edf 413.31 eps 0.0001 iteration  17
AICc 84747.15 logPost -2361226 logLik -41917.8 edf 413.30 eps 0.0001 iteration  18
AICc 84747.12 logPost -2361226 logLik -41917.8 edf 413.29 eps 0.0000 iteration  19
AICc 84747.12 logPost -2361226 logLik -41917.8 edf 413.29 eps 0.0000 iteration  19
elapsed time:  1.95min
Starting the sampler...

|                    |   0% 11.86min
|*                   |   5% 11.02min 34.79sec
|**                  |  10% 10.43min  1.16min
|***                 |  15% 10.18min  1.80min
|****                |  20%  9.89min  2.47min
|*****               |  25%  9.56min  3.19min
|******              |  30%  9.07min  3.89min
|*******             |  35%  8.64min  4.65min
|********            |  40%  8.14min  5.42min
|*********           |  45%  7.52min  6.16min
|**********          |  50%  6.88min  6.88min
|***********         |  55%  6.25min  7.64min
|************        |  60%  5.58min  8.37min
|*************       |  65%  4.90min  9.11min
|**************      |  70%  4.21min  9.83min
|***************     |  75%  3.51min 10.54min
|****************    |  80%  2.82min 11.26min
|*****************   |  85%  2.11min 11.98min
|******************  |  90%  1.41min 12.68min
|******************* |  95% 42.28sec 13.39min
|********************| 100%  0.00sec 14.46min
Code 
summary(multivariate_geo_sow_MRF_corr_endo_General)

Call:
bamlss(formula = f_rice_wheat_yield_MRF_corr_endo_General, family = mvnchol_bamlss(k = 4), 
    data = Irrig_Rev_rice_wheat, type = "modified")
---
Family: mvnchol 
Link function: mu1 = identity, mu2 = identity, mu3 = identity, mu4 = identity, lamdiag1 = log, lamdiag2 = log, lamdiag3 = log, lamdiag4 = log, lambda12 = identity, lambda13 = identity, lambda14 = identity, lambda23 = identity, lambda24 = identity, lambda34 = identity
*---
Formula mu1:
---
sowdate_fmt_rice_day ~ 1 + rice_duration_class_long + s(gw_2017) + 
    s(onset_2017) + s(monsoon_onset_dev) + s(median_onset_82_15) + 
    s(sd_onset_82_15) + s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
                            Mean    2.5%     50%   97.5% parameters
(Intercept)              71.2122 66.4929 70.9415 76.7857      2.021
rice_duration_class_long -1.3057 -1.9799 -1.3057 -0.6731     -1.311
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(gw_2017).tau21            2.250e+01 7.582e-04 2.358e+00 1.609e+02    418.029
s(gw_2017).alpha            1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(gw_2017).edf              2.288e+00 1.000e+00 1.893e+00 5.511e+00      7.436
s(onset_2017).tau21         1.829e+00 5.257e-04 6.756e-02 1.615e+01   4063.164
s(onset_2017).alpha         1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(onset_2017).edf           1.405e+00 1.000e+00 1.080e+00 3.228e+00      8.325
s(monsoon_onset_dev).tau21  5.003e+00 1.450e-04 1.193e-02 5.084e+01     14.337
s(monsoon_onset_dev).alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(monsoon_onset_dev).edf    1.616e+00 9.989e-01 1.017e+00 4.757e+00      3.543
s(median_onset_82_15).tau21 1.865e+01 1.849e-04 1.399e+00 1.308e+02     16.965
s(median_onset_82_15).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(median_onset_82_15).edf   2.491e+00 9.993e-01 1.918e+00 5.734e+00      3.660
s(sd_onset_82_15).tau21     1.273e+00 1.532e-04 8.222e-02 9.627e+00      6.198
s(sd_onset_82_15).alpha     1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(sd_onset_82_15).edf       1.366e+00 9.989e-01 1.072e+00 2.860e+00      2.558
s(District,id='mrf1').tau21 5.473e+01 1.751e+01 4.922e+01 1.163e+02   1314.536
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='mrf1').edf   3.439e+01 3.379e+01 3.443e+01 3.472e+01     34.973
s(District,id='re2').tau21  1.581e+04 9.722e+03 1.523e+04 2.483e+04      1.091
s(District,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='re2').edf    3.599e+01 3.598e+01 3.599e+01 3.599e+01     34.775
---
Formula mu2:
---
b_grain_yield_ton_per_ha_rice ~ 1 + rice_duration_class_long + 
    s(Res_rice_sow) + s(g_q5305_irrig_times_rice) + s(nperha_rice) + 
    s(p2o5perha_rice) + s(District, bs = "mrf", xt = list(penalty = K)) + 
    s(District, bs = "re")
-
Parametric coefficients:
                            Mean    2.5%     50%   97.5% parameters
(Intercept)              1.01491 0.10146 1.01126 1.77055      0.001
rice_duration_class_long 0.10392 0.03513 0.10509 0.17119      0.100
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                       Mean      2.5%       50%     97.5%
s(Res_rice_sow).tau21             3.964e+00 5.369e-02 1.530e+00 2.156e+01
s(Res_rice_sow).alpha             1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(Res_rice_sow).edf               5.470e+00 2.274e+00 5.534e+00 8.236e+00
s(g_q5305_irrig_times_rice).tau21 2.402e+00 2.424e-01 1.566e+00 1.001e+01
s(g_q5305_irrig_times_rice).alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(g_q5305_irrig_times_rice).edf   6.190e+00 4.256e+00 6.183e+00 7.877e+00
s(nperha_rice).tau21              6.460e-01 2.801e-04 2.783e-01 3.986e+00
s(nperha_rice).alpha              1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(nperha_rice).edf                3.688e+00 1.018e+00 3.711e+00 6.729e+00
s(p2o5perha_rice).tau21           2.011e-01 7.523e-05 1.685e-02 1.572e+00
s(p2o5perha_rice).alpha           1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(p2o5perha_rice).edf             2.180e+00 1.005e+00 1.703e+00 5.307e+00
s(District,id='mrf1').tau21       1.383e-02 8.222e-04 4.399e-03 7.249e-02
s(District,id='mrf1').alpha       1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='mrf1').edf         2.332e+01 1.107e+01 2.306e+01 3.291e+01
s(District,id='re2').tau21        9.451e+00 4.572e+00 8.865e+00 1.832e+01
s(District,id='re2').alpha        1.000e+00 1.000e+00 1.000e+00 1.000e+00
s(District,id='re2').edf          3.579e+01 3.564e+01 3.580e+01 3.590e+01
                                  parameters
s(Res_rice_sow).tau21                 20.478
s(Res_rice_sow).alpha                     NA
s(Res_rice_sow).edf                    8.237
s(g_q5305_irrig_times_rice).tau21      2.267
s(g_q5305_irrig_times_rice).alpha         NA
s(g_q5305_irrig_times_rice).edf        6.578
s(nperha_rice).tau21                   0.696
s(nperha_rice).alpha                      NA
s(nperha_rice).edf                     4.704
s(p2o5perha_rice).tau21                0.184
s(p2o5perha_rice).alpha                   NA
s(p2o5perha_rice).edf                  3.238
s(District,id='mrf1').tau21            0.031
s(District,id='mrf1').alpha               NA
s(District,id='mrf1').edf             31.405
s(District,id='re2').tau21             8.942
s(District,id='re2').alpha                NA
s(District,id='re2').edf              35.801
---
Formula mu3:
---
sowdate_fmt_wheat_day ~ 1 + variety_type_NMWV + s(harvest_day_rice) + 
    s(District, bs = "mrf", xt = list(penalty = K)) + s(District, 
    bs = "re")
-
Parametric coefficients:
                     Mean    2.5%     50%   97.5% parameters
(Intercept)       115.121 112.516 114.818 118.722      3.803
variety_type_NMWV  -3.058  -3.637  -3.060  -2.484     -3.071
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 1.0000 0.9999 1.0000     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(harvest_day_rice).tau21     236.645    48.673   165.224   884.720    916.881
s(harvest_day_rice).alpha       1.000     1.000     1.000     1.000         NA
s(harvest_day_rice).edf         5.642     4.260     5.588     7.310      8.328
s(District,id='mrf1').tau21   131.576    76.275   124.873   221.170   3856.743
s(District,id='mrf1').alpha     1.000     1.000     1.000     1.000         NA
s(District,id='mrf1').edf      34.777    34.648    34.783    34.874     34.993
s(District,id='re2').tau21  49601.295 30417.972 47292.497 79650.851      1.091
s(District,id='re2').alpha      1.000     1.000     1.000     1.000         NA
s(District,id='re2').edf       35.997    35.995    35.997    35.998     34.775
---
Formula mu4:
---
l_ton_per_hectare ~ 1 + variety_type_NMWV + s(Res_wheat_sow) + 
    g_q5305_irrig_times + nperha + p2o5perha + s(District, bs = "mrf", 
    xt = list(penalty = K)) + s(District, bs = "re")
-
Parametric coefficients:
                          Mean       2.5%        50%      97.5% parameters
(Intercept)         -0.5714755 -1.3450578 -0.7865354  0.9747986     -1.426
variety_type_NMWV    0.3415470  0.2931878  0.3422559  0.3875675      0.335
g_q5305_irrig_times  0.4264346  0.4003322  0.4264261  0.4545629      0.427
nperha               0.0015378  0.0009752  0.0015281  0.0021442      0.001
p2o5perha            0.0026334  0.0014479  0.0026528  0.0038092      0.003
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(Res_wheat_sow).tau21      1.964e-02 9.249e-05 2.737e-03 1.434e-01      0.181
s(Res_wheat_sow).alpha      1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(Res_wheat_sow).edf        1.618e+00 1.013e+00 1.309e+00 3.857e+00      4.461
s(District,id='mrf1').tau21 7.045e-03 8.494e-05 1.483e-03 3.990e-02      0.013
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='mrf1').edf   2.045e+01 4.016e+00 2.125e+01 3.314e+01     31.253
s(District,id='re2').tau21  5.020e+00 3.115e-01 5.093e+00 1.078e+01      4.225
s(District,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='re2').edf    3.571e+01 3.461e+01 3.585e+01 3.593e+01     35.823
---
Formula lamdiag1:
---
lamdiag1 ~ 1
-
Parametric coefficients:
              Mean   2.5%    50%  97.5% parameters
(Intercept) -2.215 -2.238 -2.215 -2.194     -2.207
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9896 0.8964 1.0000     1
---
Formula lamdiag2:
---
lamdiag2 ~ 1
-
Parametric coefficients:
               Mean    2.5%     50%   97.5% parameters
(Intercept) 0.06858 0.04750 0.06887 0.08681      0.077
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9916 0.9304 1.0000     1
---
Formula lamdiag3:
---
lamdiag3 ~ 1
-
Parametric coefficients:
              Mean   2.5%    50%  97.5% parameters
(Intercept) -2.151 -2.174 -2.151 -2.131     -2.146
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9910 0.9318 1.0000     1
---
Formula lamdiag4:
---
lamdiag4 ~ 1
-
Parametric coefficients:
              Mean   2.5%    50%  97.5% parameters
(Intercept) 0.5159 0.4954 0.5158 0.5379      0.526
-
Acceptance probability:
        Mean   2.5%    50% 97.5%
alpha 0.9906 0.9226 1.0000     1
---
Formula lambda12:
---
lambda12 ~ 1
-
Parametric coefficients:
                 Mean      2.5%       50%     97.5% parameters
(Intercept)  0.002047 -0.025822  0.003899  0.020883     -0.004
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda13:
---
lambda13 ~ 1
-
Parametric coefficients:
                Mean     2.5%      50%    97.5% parameters
(Intercept) 0.007472 0.003715 0.007414 0.011311      0.008
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda14:
---
lambda14 ~ 1
-
Parametric coefficients:
                Mean     2.5%      50%    97.5% parameters
(Intercept) 0.009566 0.004276 0.009192 0.016723      0.011
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda23:
---
lambda23 ~ 1
-
Parametric coefficients:
                 Mean      2.5%       50%     97.5% parameters
(Intercept)  0.023629 -0.008737  0.023369  0.054572      0.027
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Formula lambda24:
---
lambda24 ~ s(Res_rice_sow) + s(Res_wheat_sow) + rice_duration_class_long + 
    s(gw_2017) + s(onset_2017) + variety_type_NMWV + s(District, 
    bs = "mrf", xt = list(penalty = K)) + s(District, bs = "re")
-
Parametric coefficients:
                             Mean     2.5%      50%    97.5% parameters
(Intercept)              -0.17388 -0.24064 -0.17294 -0.10817      0.003
rice_duration_class_long -0.02739 -0.10782 -0.02621  0.04320     -0.018
variety_type_NMWV        -0.09274 -0.16866 -0.09286 -0.02381     -0.069
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
-
Smooth terms:
                                 Mean      2.5%       50%     97.5% parameters
s(Res_rice_sow).tau21       5.873e-02 7.741e-05 1.932e-03 6.186e-01      3.413
s(Res_rice_sow).alpha       1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(Res_rice_sow).edf         1.477e+00 1.003e+00 1.079e+00 4.081e+00      6.056
s(Res_wheat_sow).tau21      3.052e-02 9.663e-05 4.552e-03 2.390e-01      0.012
s(Res_wheat_sow).alpha      1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(Res_wheat_sow).edf        1.427e+00 1.004e+00 1.179e+00 3.250e+00      1.393
s(gw_2017).tau21            2.660e-02 8.469e-05 2.532e-03 1.980e-01      0.000
s(gw_2017).alpha            1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(gw_2017).edf              1.426e+00 1.004e+00 1.125e+00 3.064e+00      0.998
s(onset_2017).tau21         3.971e-02 1.315e-04 3.828e-03 3.291e-01      0.389
s(onset_2017).alpha         1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(onset_2017).edf           1.552e+00 1.012e+00 1.285e+00 3.605e+00      3.746
s(District,id='mrf1').tau21 4.122e-03 4.681e-04 3.871e-03 9.062e-03      0.004
s(District,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='mrf1').edf   1.856e+01 6.052e+00 1.950e+01 2.500e+01     20.000
s(District,id='re2').tau21  4.163e-03 7.033e-05 1.282e-03 2.427e-02      0.022
s(District,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
s(District,id='re2').edf    7.450e+00 2.723e-01 4.289e+00 2.416e+01     23.235
---
Formula lambda34:
---
lambda34 ~ 1
-
Parametric coefficients:
                Mean     2.5%      50%    97.5% parameters
(Intercept) -0.02450 -0.06433 -0.02339  0.02065     -0.008
-
Acceptance probability:
      Mean 2.5% 50% 97.5%
alpha    1    1   1     1
---
Sampler summary:
-
DIC = 84359.96 logLik = -42066.9 pd = 226.1606
runtime = 870.34
---
Optimizer summary:
-
AICc = 84747.12 edf = 413.2929 logLik = -41917.81
logPost = -2361227 nobs = 4447 runtime = 117.27
Code 
plot(multivariate_geo_sow_MRF_corr_endo_General, model = "lambda24", term = "s(gw_2017)")

Code 
plot(multivariate_geo_sow_MRF_corr_endo_General, model = "lambda24", term = "s(onset_2017)")

Code 
plot(multivariate_geo_sow_MRF_corr_endo_General, model = "lambda24", term = "s(Res_rice_sow)")

Code 
plot(multivariate_geo_sow_MRF_corr_endo_General, model = "lambda24", term = "s(Res_wheat_sow)")

7 Economics structural system analysis

8 Visualization

Code 
## | context: server
# library(distreg.vis)


# if (interactive()) {
#   distreg.vis:: vis()
# }

# library(shinylive)
# shinylive::export(appdir = "C:/Users/MMKONDIWA/OneDrive - CIMMYT/Documents/GitHub/Rice_Wheat_System_SDS_Tool/RW_System_Visualizer", destdir = "C:/Users/MMKONDIWA/OneDrive - CIMMYT/Documents/GitHub/Rice_Wheat_System_SDS_Tool/RW_System_Visualizer/docs")
#
# usethis::use_github_action(url="https://github.com/posit-dev/r-shinylive/blob/actions-v1/examples/deploy-app.yaml")


#rsconnect::deployApp("C:/Users/MMKONDIWA/OneDrive - CIMMYT/Documents/GitHub/Rice_Wheat_System_SDS_Tool/RW_System_Visualizer")

9 Conclusion

10 References